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SECTION  1 


INTRODUCTION 


The  hydrodynamics  of  estuaries  and  tidal  inlets  are  ex¬ 
tremely  complicated,  often  requiring  both  numerical  and  phys¬ 
ical  model  simulations  to  address  certain  problems  or  poten¬ 
tial  problems.  One  of  the  most  difficult  of  these  problems 
is  that  associated  with  density  stratification  in  which  fresh 
water  overrides  the  denser,  more  saline  seawater.  The  origi¬ 
nal  (December  1980)  In-House  Laboratory  Independent  Research 
(ILIR)  proposal  was  aimed  at  addressing  this  problem  of  den¬ 
sity  stratified  flow.  A  thorough  analysis  of  the  interre¬ 
lated  phenomena  associated  with  this  type  of  flow  regime 
showed  that  a  thorough  analysis  of  the  basic  turbulent  flow 
regime  associated  with  tidally  induced  flow  systems  was  first 
necessary  before  any  significant  progress  in  density  strati¬ 
fication  could  be  made.  For  example,  the  flow-induced  turbu¬ 
lence  in  a  particular  flow  regime  can  result  in  mixing  of  the 
fresh  and  saline  water,  producing  a  well-mixed  condition.  In 
this  case,  density  stratification  does  not  exist.  The  exis¬ 
tence,  partial  existence,  or  nonexistence  of  density  strati¬ 
fication  is  therefore  directly  related  to  the  basic  param¬ 
eters  of  the  flow  regime  and  to  how  these  flow  parameters 


affect  the  level  of  turbulence  in  the  flow  system.  For  these 
reasons,  the  original  scope  of  work  was  narrowed  to  an  analy¬ 
sis  of  turbulent  flow  for  a  tidally  induced  flow  system.  The 
present  analysis  does  not  include  the  effects  of  density 
stratification  since  a  complete  understanding  of  the  hydrody¬ 
namics  of  turbulent,  unstratified,  tidally  induced  flow  is 
first  necessary  before  any  meaningful  analysis  of  density 
effects  can  be  made.  Whether  or  not  density  stratification 
is  considered,  however,  the  equations  governing  turbulent 
flow  must  still  be  numerically  approximated.  If  stratified 
flow  is  then  to  be  considered,  separate  equations  must  be 
provided  to  solve  the  convection/diffusion  equation  for 
salinity  and  an  equation  of  state  must  be  provided  to  relate 
salinity  to  density.  Therefore  it  can  be  seen  that  the  pres¬ 
ent  analysis  of  turbulent  flow  is  not  inapplicable  to  strati¬ 
fied  flow  but  merely  provides  the  first,  and  most  difficult, 
step  necessary  for  analyzing  it. 

Hydrodynamic  numerical  models  approximating  tidally  in¬ 
duced  turbulent  flow  regimes  are  commonly  utilized  for  the 
evaluation  of  the  future  effects  of  some  alteration  to  the 
existing  flow  system.  Essentially  all  large  hydrodynamic 
codes  presently  in  use  approximate  the  turbulent  fluctuation 
portion  of  the  Reynolds  equations  (Turbulent  Navier-Stokes 
equations)  through  the  use  of  empirical  relationships.  Com¬ 
mon  examples  are  those  of  the  eddy  viscosity  or  "effective" 
viscosity  type  in  which  the  turbulent  velocity  fluctuations 


(Reynolds  stresses)  are  assumed  proportional  to  the  mean 
velocity  gradient  in  a  manner  analogous  to  the  viscous 
stresses.  This  is  done  through  the  use  of  an  "effective" 
viscosity  coefficient  which  accounts  for  both  the  viscous 
stresses  and  the  turbulent  stresses.  These  relationships  do 
not  realistically  account  for  the  turbulent  flow  but  instead 
employ  certain  geometry  and  flow-dependent  coefficient  values 
such  that  the  numerical  model  results  match  the  observed  pro¬ 
totype  behavior.  Thorough  descriptions  of  models  of  this 
type  are  presented  by  Launder  and  Spaulding  (1972)  and  Rodi 
(1980) .  This  approach  is  acceptable  for  many  of  the  simpler 
flow  regimes  in  which  turbulence-driven  secondary  flows  are 
not  of  primary  interest.  If,  however,  the  flow  regime  is  a 
complex  one  (such  as  the  case  of  an  estuary)  or  the  proposed 
changes  to  be  investigated  are  substantially  different  from 
the  existing  conditions,  it  can  be  seen  that  this  empirical 
approach  will  not  adequately  represent  the  turbulent  flow 
system.  The  primary  reason  for  this  failure  is  that  the 
"effective"  viscosity  models  do  not  account  for  the  transport 
of  turbulence  through  convection  or  diffusion.  When  these 
mechanisms  are  not  accounted  for,  turbulence-driven  secondary 
flows  cannot  be  predicted  to  an  adequate  degree.  If  advances 
are  to  be  made  in  the  field  of  numerical  modeling  of  estu¬ 
arine  systems,  models  which  are  capable  of  reproducing  the 
structure  of  the  turbulent  quantities  in  addition  to  the  mean 
flow  quantities  will  have  to  be  developed  and  utilized. 


There  are  presently  three  approaches  to  the  development 
of  turbulence  modeling  which  are  potentially  applicable  to 
complex  flow  regimes.  A  brief  description  of  these  models  is 
presented  below. 

1.  The  k-e  Model.  This  approach  to  the  modeling  of  turbu¬ 
lence  represents  a  very  advanced  effective  viscosity 
model  in  which  the  effective  viscosity  is  represented  in 
terms  of  turbulent  quantities.  These  quantities,  the 
turbulent  kinetic  energy  k  and  the  isotropic  dissipation 
e,  are  determined  from  separate  transport  equations.  A 
thorough  description  of  this  model  in  addition  to  vari¬ 
ous  applications  of  it  can  be  found  in  Rodi  (1980).  An 
additional  example  of  an  application  of  this  model  to  a 
long  wave  in  an  open  channel  is  presented  by  Smith  and 
Takhar  (1979) .  This  particular  model  represents  a  sig¬ 
nificant  advance  in  earlier  turbulent  models  since  the 
effective  viscosity  is  based  on  temporarily  and  spa¬ 
tially  varying  turbulent  quantities;  however,  it  does 
not  deal  with  the  actual  governing  equations  of  turbu¬ 
lent  flow.  For  this  reason,  the  applicability  may  be 
limited  to  only  the  more  basic  flow  systems. 

2.  The  Reynolds  Stress  Transport  Model.  This  approach 
attempts  to  solve,  or  approximate,  the  exact  equations 
derived  for  the  transport  of  the  individual  Reynolds 
stresses.  It  represents  the  most  rigorous  approach  that 
can  be  taken  and  is  presently  the  state-of-the-art 


numerical  treatment  of  turbulent  flow.  Because  of  the 
additional  equations  and  variables  involved,  it  also 
represents  the  most  expensive  approach.  For  this  rea¬ 
son,  the  Reynolds  stress  transport  models  have  been  ap¬ 
plied  only  to  steady-state  flow  conditions.  Excellent 
examples  of  this  formulation  have  been  presented  by 
Hanjalic  and  Launder  (1972) ,  Launder,  Reece,  and  Rodi 
(1975) ,  and  Legner  and  Finson  (1980) .  This  class  of 
models  clearly  represents  the  approach  with  the  greatest 
potential  for  modeling  flow  systems  in  which  the  trans¬ 
port  of  turbulent  quantities  is  important. 

Algebraic  Models.  A  disadvantage  of  the  Reynolds  stress 
transport  approach  is  that  multiple  partial  differential 
equations  must  be  approximated.  Some  of  these  difficul¬ 
ties  have  been  overcome  by  simplifying  the  exact  trans¬ 
port  equations  to  algebraic  relationships.  These 
approximations  generally  contain  both  the  turbulent  ki¬ 
netic  energy  k  and  the  isotropic  dissipation  e  as  unknowns 
in  addition  to  the  Reynolds  stresses.  This  approach 
therefore  represents  a  level  of  complexity  between  the 
k-e  model  and  the  Reynolds  stress  transport  model.  A 
recent  application  of  this  approach  is  presented  by  Noat 
and  Rodi  (1982)  for  steady-state  channel  flow.  Although 
this  type  of  model  is  more  economical  than  the  Reynolds 
stress  transport  models,  it  still  lacks  the  potential 
universality  of  the  transport  models. 


The  above  comparisons  strongly  indicate  that  the  best 
approach  for  effectively  dealing  with  flow  systems  as  compli¬ 
cated  as  estuaries  is  the  Reynolds  stress  transport  model. 

The  goal  of  this  project,  therefore,  is  to  develop  a  two- 
dimensional  (horizontal  and  vertical)  computational  model 
based  on  the  complete  turbulence  equations  that  will  address 
the  specific  application  of  tidally  driven  shallow-water  wave 
flow  regimes. 

A  theoretical  examination  of  the  governing  equations,  as 
they  apply  to  the  case  of  shallow-water  wave  propagation,  is 
made  to  determine  the  dominant  terms.  This  phase  of  develop¬ 
ment  is  based  on  the  concept  of  scaling  and  nondimensionali- 
zation  of  the  complete  two-dimensional  set  of  equations  for 
turbulent  flow.  The  rationale  of  this  procedure  is  begun  in 
Section  2  with  a  description  of  the  parameters  used  to  define 
shallow-water  waves.  This  definition  is  presented  to  define 
the  flow  system  for  which  the  turbulent  model  is  formulated. 
The  mechanics  of  the  scaling  and  nondimensionalization  are 
presented  in  Sections  3  and  4.  The  end  result  of  this  proce¬ 
dure  is  an  order  of  magnitude  ranking  of  the  individual  terms 
in  the  governing  equations.  This  is  done  by  first  defining 
the  relative  magnitude  of  the  individual  terms  according  to 
an  exponential  power  of  some  small  parameter  e.  This  initial 
grouping  of  terms  of  similar  magnitude  is  further  expanded 
based  on  the  assumption  that  each  flow  variable  can  be  ex¬ 
pressed  in  terms  of  a  power  series  expansion  of  the  same 


small  parameter  e.  This  formal  perturbation  expansion  tech¬ 
nique  demonstrates  the  relative  magnitudes  of  the  individual 
terms  and  provides  valuable  insight  into  the  processes  of 
turbulence.  This  insight  is  crucial  for  the  following  rea¬ 
sons:  the  equations  governing  turbulent  flow  contain  more 

unknown  quantities  than  equations.  This  indeterminant  situa¬ 
tion  can  be  referred  to  as  an  open  set  of  equations.  In 
order  to  close  the  set,  approximations  for  certain  terms  in 
the  Reynolds  stress  transport  equations  must  be  made.  These 
approximations  are  used  so  that  the  unknown  terms,  represent¬ 
ing  diffusion,  pressure  redistribution,  and  dissipation  of 
turbulent  energy,  can  be  represented  in  terms  of  the  known 
Reynolds  stresses.  The  understanding  of  the  relative  magni¬ 
tudes  of  these  terms  gained  from  the  perturbation  expansion 
analysis  is  therefore  invaluable  in  the  selection  of  an  ap¬ 
propriate  closure  approximation.  The  selected  closure  models 
and  the  derivations  of  these  models  are  presented  in  Sec¬ 
tion  5.  The  final  closed  set  of  equations  representing  the 
two-dimensional,  shallow-water  wave,  turbulent  flow  system 
are  presented  in  Section  6. 

The  application  phase  of  this  project  is  to  develop  a 
numerical  model  of  the  governing  equations  and  to  apply  it  to 
a  flow  case  for  which  complete  experimental  data  were  avail¬ 
able.  A  finite  difference  model  was  developed  for  this  pur¬ 
pose.  The  numerical  procedure  used  in  this  model  and  the 
associated  boundary  conditions  necessary  are  presented  in 


Section  7.  Verification  of  the  model  was  made  to  the  steady- 
state  flume  data  collected  by  McQuivey  in  1973.  A  descrip¬ 
tion  of  this  data  and  the  results  of  the  verification  of  the 
model  to  that  data  are  presented  in  Section  8.  In  addition 
to  the  verification  plots  and  comparisons,  an  analysis  is 
made  of  the  separate  closure  schemes  used  in  the  model. 
Finally,  Section  9  discusses  conclusions  and  future  refine¬ 
ments  and  applications  of  the  turbulence  model. 


SECTION  2 


SHALLOW-WATER  WAVES 

The  development  of  a  turbulence  model  for  shallow- 
water  waves  must  be  preceded  by  a  description  of  the 
specific  characteristics  of  those  waves.  It  must  also  be 
shown  that  these  characteristics  are  similar  to  those 
occurring  in  a  tidally  induced  environment  such  as  that 
found  in  a  typical  estuary.  This  definition  is  necessary 
for  the  development  of  nondimensionalization  parameters 
that  will  be  used  extensively  in  the  following  sections. 

In  general,  shallow-water  waves  (also  referred  to  as 
long  waves)  represent  that  class  of  waves  which  have  a 
wavelength  much  greater  than  the  depth  of  fluid  in  which 
they  are  propagating.  The  implication  of  this  statement  is 
that  the  waves  have  a  gradual  slope.  Tides,  or  astronomi¬ 
cally  induced  waves,  in  coastal  areas  generally  fall  into 
this  category  since  the  depths  of  flow  are  usually  small. 
Interest  in  the  propagation  of  small  amplitude  waves  was 
begun  in  the  nineteenth  century.  Airy  (1845)  used  the 
following  parameter  to  characterize  a  long  wave: 


where  h  is  the  depth  of  fluid  and  L  is  the  distance  between 
two  consecutive  nodes  (wavelength) .  Using  this  definition 
and  the  assumption  of  hydrostatic  pressure,  Airy  (1845) 
derived  a  set  of  equations  which  show  that  a  wave  cannot 
be  propagated  without  a  change  in  shape.  It  was  later 
shown  by  Boussinesq  (1871)  and  independently  by  Rayleigh 
(1876)  that  a  permanent  solitary  long  wave  did,  in  fact, 
exist  if  the  effects  of  dynamic  pressures  were  accounted 
for.  This  contradition  of  whether  or  not  a  long  wave  can 
propagate  without  change  of  shape  led  to  the  "long-wave 
paradox"  disagreement  reported  by  Lamb  (1932)  and  later  by 
Ursell  (1953)  and  Madsen  and  Mei  (1969) . 

Friedrichs  (1948)  presented  a  derivation  of  the 
governing  equations  which,  while  not  giving  a  complete 
explanation  of  the  paradox,  began  to  shed  considerable 
light  on  the  problem.  Friedrichs'  (1948)  derivation  was 
based  on  a  perturbation  expansion  in  which  the  expansion 
parameter  e  was  similar  to  that  defined  by  equation  2.1. 

He  showed  that  the  first  approximation  to  the  equations 
was  consistent  with  the  results  of  Airy,  i.e.  a  wave  of 
permanent  form  did  not  exist.  A  second  order  approximation 
of  Friedrichs'  derivation  by  Keller  (1948),  however,  did 
show  the  existence  of  a  wave  of  permanent  shape.  Even 
higher  order  approximations  of  the  permanent  wave  were 
presented  by  Wehausen  and  Laitone  (1960)  by  extending 
Friedrichs'  approach. 


The  first  satisfactory  explanation  to  the  long-wave 
paradox  was  presented  by  Ursell  (1953)  .  He  demonstrated 
that  the  parameter 


where  ^  is  a  measure  of  wave  amplitude,  should  be  taken 
into  consideration  in  the  derivation  of  long-wave  theory 
as  a  "fundamental  critical  parameter."  His  derivations 
showed  the  following  three  cases  indicating  three  types  of 
long  waves 

Case  Is  e<<l.  The  governing  equation  becomes  linear 
with  no  permanent  form  solutions. 

Case  2:  e=0(l).  The  governing  equation  were  of  the 

Boussinesq  type  and  therefore  showed  permanent 
form  solutions. 

Case  3:  e>>l.  The  governing  equations  become  of  the 

Airy  type  with  no  permanent  form  solution. 

The  above  descriptions  were  presented  so  that  the 
proper  scaling  parameters  can  be  utilized  for  the  non- 
dimensionalization  of  the  governing  equations  for  the  appli¬ 
cation  to  estuarine  conditions.  At  this  point,  typical 
values  for  an  estuary  should  be  defined.  The  following 
values  are  approximate  averages  for  the  Chesapeake  Bay  but 
are  representative  of  other  estuaries. 

C  «  2  ft 
L  200  miles 


These  values  yield  the  following  values  for  the  expansion 
parameters  shown  in  equations  2.1  and  2.2 


and 


The  above  values  demonstrate  that  typical  estuarine  systems 
are  governed  by  Airy  type  equations.  Consideration  of  this 
result  led  to  the  conclusion  that  the  Friedrichs'  expansion 
parameter  should  be  a  valid  one  for  use  in  a  tidally  induced 
estuarine  system.  As  previously  mentioned,  Wehausen  and 
Laitone  (1960)  used  a  similar  expansion  parameter  in  their 
extensive  expansion  of  the  shallow-water  wave  equations 
for  the  case  of  inviscid,  irrotational  flow.  Since  this 
approach  has  shown  to  be  valid  for  that  case,  it  should  be 
equally  valid  for  the  case  of  turbulent  shallow-water  wave 
propagation.  The  following  two  sections  will  therefore 
make  use  of  Friedrichs'  expansion  parameter  to  develop  the 
governing  equations  for  viscous,  rotational,  turbulent 


£  *  2.8  x  10~5 
Li 


«  8. 3  x  10' 
.3 


shallow-water  waves. 


SECTION  3 


THE  REYNOLDS  EQUATIONS  AND  THE  CONTINUITY  EQUATION 

The  governing  equations  for  a  constant  density,  incom¬ 
pressible,  Newtonian  fluid  are  the  Navier-Stokes  equations 
and  the  continuity  equation.  The  turbulent  equivalent  of 
these  relationships  are  the  Reynolds  equations  and  the 
turbulent  continuity  equation.  The  two-dimensional  form 
can  be  written  as  follows: 


Reynolds  Equations 


9U 

9t 


u  i  +  v 


9V 

9y 


1  9P 
p  9x 


2— 
9  U 

9  2 
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9u 1 u 1  -  9u' v 
9x  9y 


3.1 


9X  +  U9V  9V  =  _19P 

9t  9x  v  9y  p  9y 


glTy1-  _  dvTvT 
9x  9y 

3.2 


Continuity  Equation 

9U  9V  _  _  3.3 

9x  9y  u 


In  addition  to  the  above  relationships,  the  boundary  condi 
tions  for  the  free  surface  and  bottom  can  be  written  as: 


Surface  Boundary  Equation 


In  the  above  equations,  the  overbars  represent  the  averaged 
value  of  either  the  mean  flow  quantity  (in  upper  case)  or 
the  turbulent  fluctuation  (represented  by  the  primed  quan¬ 
tities)  .  The  individual  terms  are  defined  as 
U  -  mean  horizontal  velocity 
V  -  mean  vertical  velocity 
P  -  mean  pressure 

u 'u* ,  u' v' ,  v'v'  -  Reynolds  stresses  (mean 

value  of  turbulent  velocity  fluctuation 
correlation) 

^  -  water-surface  elevation  above  some  datum 
H  -  depth  of  flow 
x  -  primary  horizontal  axis 
y  -  primary  vertical  axis 
v  -  kinematic  viscosity 

A  complete  derivation  of  equations  3.1  through  3.5  is  pre¬ 
sented  in  appendices  A  and  B. 

In  order  to  begin  the  order  of  magnitude  analysis,  it 
is  necessary  to  nondimens ionalize  the  governing  equations 
such  that  the  horizontal  and  vertical  distances  are  stretched 


by  different  amounts.  The  following  nondimensional  terms 
are  therefore  defined: 
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3.6 


3.7 

3.8 

3.9 

3.10 

3.11 

3.12 

3.13 


The  above  definitions  make  use  of  the  Friedrichs'  expansion 
parameter  such  that 


£ 


1 


where  h  is  a  length  scale  representing  depth  of  flow  and  1 
represents  a  wavelength  scale.  The  reference  velocity  VQ 
represents  the  velocity  of  a  gravity  wave  in  a  fluid  of 
depth  h.  The  turbulent  fluctuation  terms  appearing  in 


equations  3.1,  3.2,  and  3.4  are  nondimensionalized  by  use 
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of  an  as  yet  undetermined  coefficient.  These  relationships 
are  written  as  follows: 


u;  Hi 

9x 


u ' 


Hj 

9x 


/Y 


3.14 


and 


u  f  u !  =  u  !  u /a 
i  j  13  x 


I 

3.15  I 


where  indicial  notation  has  been  used  for  the  Reynolds 
stresses  u ' u 1 ,  u ' v ' ,  and  v 1 v ' .  A  reference  Reynolds  number 
can  also  be  defined  based  on  the  reference  velocity. 


R  = 


Voh 


3.16 


Substitution  of  the  nondimensional  relationships  of  3.6  - 
3.15  into  equations  3.1  -  3.5  and  grouping  according  to 
powers  of  the  expansion  coefficient  e  =  h/1  yields  the 
following  equations: 


0  ~9U  1  :2U  .  / 9U  p9U  3P, 

£  (V— 7T  -  p  >  A  J  )  +  £  {—7T  +  U~ +  —>r) 

9y  dy^  9 t  3x  9x 
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+  e  ( — ~  — - — )  +  e  ( — ^  — * — ) 

V  3x  V  9y 

o  o 


e2  (± 

'R 

9x 


=  0 


3.17 
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A  A  a  A  A  A  a  A 

0  ,r.3V  3P  1  sV  .  ,3V  r.3V.  2,1  3  v. 

£  (V— X-  +  -7T  +  1  -  p  —>r~)  +  £  (— X-  +  U — x-)  -  £  (p  — x^) 

3y  3y  *  Dy ^  3t  3x  K  3xz 


,  ,  a  3u'v'.  ,  o.  a  3v'v'. 

+  £  ( — ~  — x — )  +  £  ( — j  — x — )  =  0 

V  z  3x  V  z  dy 
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3.18 


,3U,  .  o,3V.  . 

£  (— )  +  e  (— xr )  =  0 
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£  (—5-  +  u  —%•)  +  £°  u^5l  -  V  > 
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3.20 


£  (—S’)  +  V  =  0 
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3.21 


An  analysis  of  this  preliminary  grouping  leads  to  some 
simplifications  which  can  be  used  in  the  further  expansions 
of  the  terms.  For  example,  the  following  conclusions  can 
be  deduced: 

a)  equation  3.19  shows  that 


0(e) 


3y 


or  that  the  vertical  velocity  V  is  essentially  constant 


throughout  the  water  column.  Equation  3.21,  however,  shows 
that  the  bottom  velocity  is  small;  for  example: 


It  can  then  be  seen  that  the  constant  vertical  velocity  is 
small  and  can  be  expressed  as  shown  below 


V  =  o(e) 


b)  The  as  yet  undefined  coefficient  a  must  be  selected 
such  that  the  final  set  of  equations  will  accurately  re¬ 
present  the  complete  range  (manifold)  of  solutions  to  the 
given  problem  of  small  amplitude  waves.  An  extreme  case 
can  be  presented  as  steady  uniform  flow  in  the  x  -  direction 
For  this  case,  equation  3.17  can  be  written  as 
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Two  conditions  can  now  be  considered.  First,  if  the  verti- 

/v 

cal  velocity  V  is  assumed  zero  and  the  flow  is  considered 

/v 

nonturbulent  (u’v1  =  0) ,  then  it  can  be  seen  that  the 
above  relationship  should  be  written  as 


2" 

3  U  3P. 

~~*~2  +  ~~K‘ 
3y  3x 


0 


which  is  the  equation  for  laminar  flow  between  parallel 
plates  (Daily  and  Harleman,  1966,  p.116).  If,  on  the 
other  hand,  turbulence  is  considered,  it  can  be  seen  that 
the  following  relationship  must  be  valid 
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3uV) 

ay 


0  (c) 
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This  leads  to  the  definition  of  the  Reynolds  stress  co¬ 
efficient  as 


3. 


ct  =  eV 

o 

Equations  3.17  through  3.19  can  now  be  written  entirely 
in  terms  of  £  as  shown  below. 
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The  above  equations  now  represent  the  initial  non- 
dimensionalization  of  the  governing  equations.  The  bound 
ary  condition  equations  3.20  and  3.21  will  be  considered 
at  a  later  point  since  they  do  not  affect  the  further 
development  of  the  Reynolds  equations  or  the  continuity 
equation. 


In  order  to  complete  the  order  of  magnitude  analysis 
of  equations  3.1  -  3.3,  it  is  necessary  to  assume  that 
the  variable  quantities  can  be  represented  in  terms  of  a 
power  series  expansion  with  respect  to  some  small  param¬ 
eter,  in  this  case  the  expansion  parameter  e  (Stoker, 

1957,  and  Ippen,  1966).  This  identical  procedure  was  used 
by  Wehausen  and  Laitone  (1960)  in  their  derivation  of  the 
nonturbulent ,  inviscid  shallow-water  wave.  Following  the 
above  examples,  the  variable  expansions  take  the  follow¬ 
ing  form: 

A  A 

U  =  V  +  eU,  +  e  u~  +  .  .  .  emU 
o  1  2  m 

V  =  eV  +  e2V2  +  .  .  .  cmVm  (since  V  =  0  (e)  3.26 

P  =  P  +  eP.  +  e2P_  +  .  .  .  emP 
o  1  2  m 

u'u'.  =  (u'u!)  +  e  (u'u'.).,  +  £2(u'u!)_  +  .  .  .  em  (u'u'. ) 

~  D  DO  -  D  1  -  D  2  ~  D  m 

Substitution  of  equations  3.26  into  equations  3.23  through 
3.25  and  collection  according  to  powers  of  e  give  the 
following  order  of  magnitude  relationships  (for  brevity, 

3 

only  the  relationships  through  0(s  )  will  be  written.  The 
overmarks  indicating  nondimensional  x  and  y  directions 
will  also  be  omitted) : 


3V,  3U. 
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3y  l3x 


3V-,  -  3U_ 

+  r-^l  +  e  [-5-=- 

3y  1  dx 


+  0  (e  )  =  0 


3.29 


Inspection  of  the  order  of  magnitude  relationships 
presented  in  equations  3.27  through  3.29  demonstrate  the 
following : 

a)  The  zero  order  solution  of  the  Reynolds 
equations  is  that  of  hydrostatic  pressure. 
Although  this  is  a  common  assumption,  the 
method  of  analysis  does  confirm  the  rela¬ 
tionship. 

b)  The  Reynolds  stress  terms  are  of  the  same 
order  of  magnitude  as  the  terms  represent¬ 
ing  convection,  pressure,  viscosity,  etc. 

The  proper  modeling  of  a  turbulent  flow 
system  then  requires  the  proper  modeling 
of  the  Reynolds  stresses. 

c)  The  system  of  equations  does  not  represent 
a  closed  set  of  equations;  they  now  contain 

the  unknown  Reynolds  stresses  u'u!. 

~  3 

Conclusions  a)  and  b)  are  evident  from  the  presented 
analysis;  however,  the  approach  to  conclusion  c)  represents 
the  difference  between  proper  turbulence  modeling  and 
"pseudo"  turbulence  modeling.  The  common  procedure  at  this 
point  is  to  represent  the  Reynolds  stresses  through  some 

22 


simplified  closing  relationship  using  parameters  such  as 
mixing  lengths  or  combinations  of  mixing  length  and  turbu¬ 
lent  kinetic  energy.  Early  models  of  this  type  have  been 
developed  by  Prandtl  (1925,1945)  and  Kolmogorov  (1942). 

Many  of  the  more  recent  closure  techniques  are  presented  by 
Launder  and  Spalding  (1972)  and  by  Rodi  (1980) .  The  ap¬ 
proach  selected  for  this  analysis  is  the  approximation  of 
the  Reynolds  stresses  based  on  the  exact  equations  for  the 
Reynolds  stresses  (Reynolds  stress  transport  equations) . 
Although  these  equations  also  do  not  represent  a  closed 
set  of  equations,  closure  at  this  level  represents  a  higher 
order  of  accuracy  turbulence  model  known  as  a  second-order 
closure  model.  The  closure  approximations  employed  for 
the  model  are  presented  in  Section  5  following  a  detailed 
analysis  of  the  Reynolds  stress  transport  equations  in  the 
following  section.  This  is  done  in  an  analagous  manner  to 
that  in  which  the  perturbated  Reynolds  equations  of  this 
section  were  developed. 


SECTION  4 


THE  REYNOLDS  STRESS  TRANSPORT  EQUATIONS 


The  governing  equations  for  the  transport  of  the  turbu¬ 
lent  Reynolds  stresses  for  a  constant  density,  incompressi¬ 
ble,  Newtonial  fluid  can  be  written  in  indicial  notation  as 
follows: 
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4.1 


The  overbarred  quantities  represent  the  averaged  values  of 
either  the  mean  flow  variables  (in  upper  case  letters)  or 
the  fluctuating  variables  (represented  by  the  primes) .  The 
above  equations  are  derived  through  the  manipulation  of  the 
Reynolds  equations  and  the  continuity  equation.  A  complete 
derivation  of  equation  4.1  is  shown  in  Appendix  C. 

Equation  4.1  represents  six  equations  for  the  Reynolds 
stresses  u ' u ' ,  v ' v ' ,  w'w' ,  u ' v ' ,  u'w' ,  and  v' w' .  The  case 

flow  (for  which  =  f  (x£)  ,  ^  =  Uj  *  0 


of  two-dimensional 


and  symmetry  of  flow  with  respect  to  planes  perpendicular 


to  the  x^  axis)  reduces  the  number  of  equations  to  four 
(uT"ur,  v' v'  ,  w'w' ,  and  u '  v ’ )  .  This  is  due  to  the  fact  that 
all  correlations  involving  u^  and  all  uneven  derivatives 
with  respect  to  the  x3  direction  are  zero  (Hinze,  1959, 
p.252,  and  Cebeci  and  Smith,  1974,  p.83).  These  four  equa¬ 
tions  now  represent  the  governing  equations  for  the  trans¬ 
port  of  the  Reynolds  stresses  for  the  two-dimensional  flow 
case.  Each  equation  can  now  be  nondimensionalized  in  a 
similar  manner  to  that  performed  for  the  Reynolds  equations 
in  the  previous  section. 

The  nondimensionalization  process  begins  with  the  use 
of  equations  3.6  through  3.10.  The  z  -  axis  direction,  now 
represented  in  equation  4.1,  is  scaled  identically  to  the 
x  -  axis  as  shown  below 

z  =  z/1  4.2 

Use  will  also  be  made  of  the  results  shown  by  equation  3.22 
for  the  Reynolds  stress  scale.  This  was  shown  to  be 

A 

u  !  u  !  =  uTuT  /a 
13  13' 


or 


uiu j  =  Avo‘ 


In  addition  to  the  above  parameters,  it  is  also  necessary 
to  develop  scaling  relationships  for  the  new  correlation 
quantities  in  equation  4.1  which  represent  pressure 


redistribution,  dissipation,  and  diffusion  of  turbulence. 
Since  each  of  these  terms  are  of  unknown  relative  magnitude 
they  will  be  addressed  separately. 


a)  Diffusion.  Scaling  for  the  triple  correlation 
term  representing  diffusion  of  turbulence  is 
based  on  the  assumption  that  a  single  fluctua¬ 
tion  can  be  written  as  follows: 

*.  _  ui 


4.3 


This  represents  an  extension  of  equation  3.22. 
Use  of  this  relationship  yields  the  following 
nondimensional  form  for  the  diffusion  term 
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4.4 


b)  Pressure  Redistribution.  The  terms  in  equation 
4.1  which  contain  pressure  fluctuation  correla¬ 
tions  can  be  consolidated  as  follows: 
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An  as  yet  undetermined  coefficient  can  now  be 
introduced  as 


U!  _Rl  =  IU  lEl  /  — _ 

j  x^  p  j  3x^  '  x^  -  scale 

where  the  term  x^  -  scale  refers  to  the  appro¬ 
priate  x,  y,  or  z  direction  scale  and  use  has 


been  made  of  equation  4.3.  A  Poisson  equation 
for  the  single  pressure  fluctuation  can  be  written 
as  follows 
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The  derivation  of  this  relationship  is  based  on 
the  divergence  of  the  Navier-Stokes  equations 
and  is  presented  in  Appendix  D.  If  the  spacial 
derivatives  are  neglected,  it  can  be  seen  that 
the  lowest  order  of  magnitude  representation  of 
the  pressure  fluctuation  can  be  written  as  follows 
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The  term  of  interest  must  then  be  expressible  as 
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This  assumption  leads  to  the  assumption  that 

0  = 


or  that  the  scaled  term  can  be  written  as 
~  _  £V  ^ 

u.  1L1  =  1  u.  lEl  /  _2 _ 

j  3x^  p  j  3x^  xi  “  scale 


A  comment  can  now  be  made  concerning  the  term 
labeled  x.  -  scale.  The  pressure  fluctuation 


term  represents  the  sum  of  the  following  two  parts 


u 


Du'p 

~Jx~ 


1 


I 


-  p' 


3u 
_1 
ax . 


It  is  generally  acknowledged  that  the  first  term 
on  the  right  is  considerably  smaller  in  magni¬ 
tude  than  the  second  term  (Hanjalic  and  Launder, 
1972) .  In  this  case,  the  scaling  must  be  based 
on  careful  examination  of  that  term.  For  example, 
in  the  equation  for  u'u',  v'v1,  and  w'w',  the 
second  term  becomes 


3x  ' 


,  3v 
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or 


which  are  all  interrelated  due  to  the  continuity 
equation.  This  would  appear  reasonable  since 
pressure  exerts  in  all  directions  equally.  With 
this  analogy  in  mind,  the  x^  -  scale  was  selected 
as  that  of  the  y  -  direction.  Therefore,  the 
final  scaling  for  the  pressure  redistribution 
term  is  as  shown  below 


A 


4.5 


c)  Dissipation.  The  terms  in  equation  4.1  which 
represent  dissipation  of  turbulence  can  initially 
be  determined  following  the  expansion  of  equation 


4.1  based  on  the  already  determined  coefficients 
This  temporary  relationship  is  as  follows: 


3x^  k  k  (x^  -  scale) 


The  individual  equations  representing  the  transport  of 
the  Reynolds  stresses  can  now  be  written  and  nondimension- 
alized  by  use  of  equations  3.6  -  3.10,  4.2,  and  4.4  -  4.6. 
The  dimensional  form  (presented  for  reference  purposes) 
and  the  nondimensional  form  are  shown  as  follows: 
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3x  3y  p  3y  9x2  gy2 


+  2v^^+2v-^-^-+2v-^-^-=0 
V  3x  3x  3y  3y  3z  3z  u 


v'v'  : 


e  +  2  uru.  IV  +  5  8VW1,  +  2v?v.  3V  +  XV. 

3t  3x  3x  3y  3y 


+  e3/2  (9u-:-v-^)  +  e1/2  L  +  2  V  *$L 

3x  3y  3y 


,  v  %  32v'v'  ,  v  ,  32v'v'  ,  _  ,2vg,  3v'  3v' 

_  £  (y  1  /  *  5  _  l \7  U  )  A  n  ^  Wr  1  *  A  ~ 

vo‘L  3xz  von  3y  vox  3x  3x 


,2v6  3v '  3v '  ,2v6,  3v'  3v' 

*  'V  h '  ^  ^  ^  ^  V  1  A  ^ 

on  3y  3y  c>  3z  3z 


4.10 


1 vmx9&^&xr&irsc 


3x  9x 


u  '  v '  : 


c  (  IH^X  +  u'v '  ™  +  u'u'  i£  +  G  iiX)  +  v'V  is 
3t  3x  3x  3x  3y 


;  ,  av  ,  3uv  ,  _ 

+  u'v*  —  +  V  * +  £ 

3y  3y 


r.  3u 1  v 1  A  3/2  3u  1  v 1  v  1  ^  1/2  Bu'v'v1 


2  ~ 


+  v.  3£l  +  u,  iEl  _  e  (_H_)  3  u'XL  _  (_v_) 

'N  7^  lV  l1  ~2  lV  hJ 


v  %  32u'v' 


o  3x 


,  .  ,2vB,  3u'  3v'  ,  ,2vBn  3u‘  3v'  ,  ,2vBx  3u'  3v'  „ 
+  e  (y-j)  -x-.+  (y-^)  7*-  +  £  IfT1  7^  7~~  =  0 


r  ii  A  « 

o  9y 


'v  V 

ox  3x  3x 


\Tr  v_  /  ^  ^ 

o  3y  3y 


o  3z  3z 
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The  final  remaining  nondimensionalizing  coefficient 

B  can  be  determined  by  inspection  of  the  zero  order  terms 

of  equation  4.10.  The  sum  of  these  terms  must  be  at  most 
1/2 

of  the  order  £  .  This  relationship  is  shown  below 


£°  [2v’v'  ^  +  V 


3V  ,  3v'v'  1  32v'v'  ,  2B  av’^Sv’  ,  ....Sp', 

/%  '  V  n  A  '  p  A  X  '  V  V  A  J 

3y  3y  3y  3y  3y  3y 


-  0  U1/2) 


From  previous  analyses,  it  has  been  shown  that: 

a)  The  first  two  terms  are  of  order  e  since  it  has 

A 

been  determined  that  V  =  0  (£) . 

2  ~ 

1  9  1 v 1 

b)  The  term  ^  ^ —  is  of  order  e  since  it  has  been 

3y 

shown  in  Section  3  that  the  term  1/R  is  small  and 
should  be  expressed  in  the  form  e  (-4)  for  proper 
scaling . 
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c)  The  pressure  fluctuation  terms  are  approximately  of 

~  OTT 

the  order  of  magnitude  of  u!u'.  -r —  and  are  therefore 

1  D  3xi 

not  insignificant  in  value. 

In  view  of  the  above  observations,  it  is  reasonable  to  as¬ 
sume  that  the  last  two  terms  are  of  equivalent  orders  of 
magnitude  for  the  following  reason:  The  pressure  fluctua¬ 
tion  term  represents  a  source  of  turbulent  energy  trans¬ 
ferred  into  the  system.  This  source  of  energy  must  be 
balanced  by  the  dissipation  term  since  the  remaining  terms 
have  been  shown  to  be  small.  The  following  relationship  is 
then  based  on  this  assumption 

v."i£l  *  u;v.  4  =  1 

3y  3y  R  3y  3y 

In  order  for  this  relationship  to  occur,  the  term  6/R  must  be 
of  order  unity.  Therefore,  the  final  scaling  parameter  is 

•-* 

Equations  4.8,  4.10,  4.12,  and  4.14  can  now  be  written  in 
a  manner  compatible  with  equations  3.23  through  3.25.  The 
parameter  1/eR  is  used  in  the  appropriate  locations  similar  to 
the  earlier  derivations.  The  following  relationships  result 


e°(2u' 


v.  ™  +  v  -3-u-Pl  +  i!il  ISil  +  2u*  +  e1/2  (lu-lu-'vl 

eR  3y  3y 


+  e  +  2u ' 


,  3U  *  3u '  u  1  1_  32u'u' .  3/2  Su'u'u' 

U  »+u  -  fR  .9  '  +  E 

3x  3x  3y  3x 


2  ,2  3u 1  3  u  '  .  2  3u '  3u' 


1  32u'u' 


v '  v '  : 


av  (  ;  av'v' 

+ 
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av ' 

ay  ay 

eR 

ay 

<>i 

cx> 

+  E1/2  until  +  e , 3v£vj_  +  2u;v,  av  +  -  9V^1  .  1  Av, 

a  ..  i.  r\  _  _  r.  __  CR  n 


ay 

at 

ax 

ax 

3U* V' V* 

,  2  ,  2 
+  E  {7r 
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,  2  av' 
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A 

A 

+  eR 

* 

ax 

ax 

ax 

az 

az 

3  ,  l  a2v'v\  _  „ 

+  E  '•  7r  — “ 1  *  0 

3x 
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+  e  ( 


3w'w'  2 

A  +  eR 

3w ' 
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+  2w ' 

i£l)  + 

ay 

ay 

ay 

3z 

3w'w'  +  n  ■ 

aw' w' 

■)  + 

3/2 

j au'w'w' 
ax 

-  i  VJ 

at 

ax 

e 

2  f—2  5wf  3w*  __2  3w»  aw'  _  n 

^  gR  a  *  eR  *  ~ 
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A 

r'v' 

3U  +  u'v'  3V  +  v  3u'V’ 
ay  ay  ay 

2  au' 

ER  ay 

av' 

ay 

+  u' 

l£l 

ay 

i£l, 

3x 

+  ei/2  +  £ 

ay 

j  a  u '  v ' 
at 

ii  1 1?  i 

au 

U  V 

3x 

1Y  + 

u  3u'v' 

1 
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3x 

8x 

eR 

ay 

2  au' 
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2 

3u' 

av' 

cR  ax 

~  T 

A 

ax 

eR 

az 

az 
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,  ,  3,  1  32u'v\  _  „ 
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In  a  manner  consistent  with  the  expansion  of  the  Reynolds 
equations,  it  is  assumed  that  the  flow  variables  can  be 
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represented  as  a  series  expansion  of  the  parameter  e. 
These  expansions  take  the  following  form 


U  =  U  +  eU,  +  £2U„  +  .  .  .  +  emu 
ol2  m 

V  =  eV1  +  £ 2V 2  +  ...  +  e1^  (since  V  =  0  (e) 


+  e1//2  (u[u;!)1  +e(u[u^)2  + 


+ 


£m/2 


(u'.u! ) 

l  j' 


m 


uiujuk  =  (uiujuk>o  +  e1/2  luiujuk>l  +  •  •  • 


+  em/2  (u!u! u') 
i  j  km 


3u^  3u ^  3u^  3Uj 

- A  ~~  ~  ^ 3x,  3x,  ^o  + 

3xk  3xk  k  k 


£l/2  + 

(3xk  3xk>l 


3ui  aui 


/0  3u .  3u . 

,  m/  2  ,  i  j . 

+  £  ^  ^,ra 


u!  Iftl- 

1  3x . 

1 


.  ,  3p 1 .  ,  1/2  ,  ,  3p  , 

(u!  Tf— )  +  £  '  (u!  -rf— )  ,  + 

i  3x .  o  i  3x .  1 


+  £m/2  (u!  |El) 
i  3x  .  m 
J 


The  expansions  of  the  fluctuating  quantities  of  equations 

1/2 

4.19  have  been  made  in  terms  of  £  This  was  done  to 
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reflect  the  fact  that  equations  4.15  through  4.18  were 

derived  by  multiplication  of  the  Reynolds  equations  (which 

were  expanded  in  terms  of  e)  by  a  velocity  fluctuation 

1/2 

u^  (which  is  of  order  c  ) .  This  expansion  should  there¬ 
fore  achieve  compatibility,  with  regard  to  magnitudes, 
with  the  expanded  Reynolds  equations.  Substitution  of 
equations  4.19  into  equations  4.15  through  4.18  and  collec¬ 
tion  according  to  powers  of  e  give  the  following  order  of 
magnitude  relationships  (for  the  sake  of  brevity,  only 
those  terms  up  to  and  including  0  (c)  will  be  shown.  The 
overmarks  indicating  nondimensional  x  and  y  directions 
will  also  be  omitted) : 


o  r~,  ,  o,2  ,9u'  9u'.  ,  ,9p' 

£  [2<u  v  >o  3?-  +  a?  (3^  Sy^o  +  2|U 


1/2  ,,  3Vo  ,  2  .  9u '  30'.  ,  S(u'u'v,)o 

+  E  [2(u  v  )2  IT  IT-1!  +  — 55 - 


a_,  au  au, 

+  2(u-  |£-);l]  +  e  [  2  (u  ’  v 1  )  2  +  2(u'V)0  ^ 


a(u’u’)  2  3  .  3  ,  gp, 

+  vi  — ay -  +  7r  {W~  +  2(u  a x  >2 


3(u,u'v’)1  3(u'u')_  9U  a(u'u') 

*  - 5^  +  at  +  2(u'u>>o  S2  +  uo  — T5T- 


1  92(u'u')c 

ER  3y2 


]  =  0 
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u  '  v  '  : 


C  °  [  ( V 1  V  ' )  —  +  ^p(—  — )  +  u'  4^-)  +  (v'  ] 

o  dy  eR  -y  dy  o  dy  o  '  3x  ;oJ 


wo  3(J  t  ^  ,  3  (u‘  v*  v* ) 

+  c 17 2  [  ( v  ’  v 1  )  —  +  IXl)  +  _ 2 

'  1  >>y  eR  1  3y  9y  '  1  3y 


+  <u'  §£->,  +  (V  Ifi.j 


3V, 


3U 


,  .w  3U.  9(u'v') 

+  £l<uV)o  ^  +  (VV)2  ^  +  (VV)0  ^  +  V,  jj—2- 


")  3u'  3v'  ^  ^  V  V  )  .  Pin1  3n* 

eR ' 3y  3y  ;2  3y  +  (u  3y  ]2  {V  3x  ]2 


3  (u 1  v 1 ) 


3U_ 


3  (u  1  v 1  ) 


3 1  +  (u'v,)o  3x  +  Uo  3x 


1  32<“,V>0 

eR  3y2 


]  =  0 
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The  general  validity  of  the  above  order  of  magnitude 
analysis  can  be  qualitatively  determined  by  comparison  with 
published  energy  balance  plots  based  on  experimental  bound¬ 
ary  layer  data.  An  assessment  of  the  equations  can  be  made 
by  analyzing  the  zero  order  terms  of  equations  4.20  through 
4.23  as  they  apply  to  the  published  test  cases  for  two- 
dimensional  flow  (U  =  f  (y) ,  V  =  0) .  For  example,  it  can  be 
seen  that  production  terms  only  occur  in  the  equations  for 
u ' u '  and  u'v'  and  that  this  production  is  balanced  by  a 
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combination  of  dissipation  and  pressure  redistribution.  All 
diffusion  terms  can  be  seen  to  be  of  a  lower  order  of  sig¬ 
nificance.  The  remaining  two  equations  (v 1 v 1  and  w'w ' ) 
have  no  production  term  but  instead  show  that  the  energy 
entering  the  equations  through  pressure  redistribution  is 
balanced  by  the  dissipation  term.  As  in  the  other  equations, 
the  diffusion  terms  are  of  less  importance.  These  general 
observations  are  verified  by  the  energy  balance  plots  of 
Townsend  (1956)  and  Laufer  (1954)  as  presented  by  Hinze 
(1959,  p.  498,  p.  531).  These  comparisons  demonstrate  that 
the  scaling  parameters  used  were  valid  in  that  the  final 
order  of  magnitude  relationships  derived  in  this  section  are 
qualitatively  confirmed  by  experimental  data. 

An  evaluation  can  now  be  made  of  the  significance  of 
the  approximations  needed  to  represent  the  unknown  terms 
with  terms  involving  Reynolds  stresses  and  mean  velocity 
values.  Based  on  the  results  of  the  order  of  magnitude 
analysis,  it  can  be  seen  that  both  the  dissipation  terms  and 
the  pressure  redistribution  terms  are  significant  in  rela¬ 
tion  to  the  production  terms.  The  approximations  for  these 
terms  must  therefore  be  carefully  considered  if  any  effective 
turbulence  model  is  to  result.  On  the  other  hand,  the  triple 
correlation  terms  representing  diffusion  can  be  approximated 
to  a  lower  degree  of  accuracy  since  these  terms  have  been 
shown  to  be  less  vital  in  the  ultimate  solution  of  the 
governing  equations.  The  following  section  provides  a 


SECTION  5 


CLOSURE  APPROXIMATION  SCHEMES 


The  importance  of  an  adequate  and  appropriate  closure 
approximation  for  the  dissipation,  pressure  redistribution, 
and  diffusion  terms  in  the  Reynolds  transport  equations 
was  stated  in  the  previous  section.  The  closure  approxima¬ 
tion  schemes  for  this  turbulence  model  were  selected,  or 
derived,  so  that  difficult  flow  situations  could  be  properly 
represented.  The  following  subsections  provide  a  detailed 
explanation  of  the  various  closure  schemes  used  and  the 
reasons  for  their  selection. 


Dissipation  Closure  Approximation 
The  terms  representing  dissipation  of  turbulent  energy 
can  be  written  as  shown  below 


The  usual  closure  approach  is  based  on  the  assumption  that 
turbulence  is  isotropic  and  that  dissipation  is  then  also 
isotropic.  This  leads  to  the  following  approximation 


2 

3 


1  3x, 


E 
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where  is  the  Kronecker  delta,  defined  as 

5ij  =  1  for  i  =  j 
6ij  =  0  for  1  ?  j 

and  e  represents  the  isotropic  dissipation  term.  This 
approach  is  possibly  valid  for  high  Reynolds  number  flows, 
but  may  possibly  lead  to  considerable  error  for  reversing 
flow  situations  common  to  tidal  estuaries.  In  view  of  the 
importance  of  the  dissipation  terms  in  the  governing  equa¬ 
tions  and  the  belief  that  dissipation  is  not  isotropic  for 
the  specific  case  of  tidally  induced  flows,  a  different 
approach  will  be  presented.  The  basis  of  this  approach  is 
the  two-point  velocity  correlation.  A  detailed  development 
of  the  scheme  is  presented,  first  for  the  auto  correlations 
(i  =  j)  and  then  for  the  cross  correlations  (.i  7*  j). 

Auto  Correlations 

The  auto  correlation  portion  of  the  dissipation  of 
turbulent  energy  term  can  be  written  as 

3u!  3u! 

( _ i.  - i) 

3xk  SV 

Since  the  following  development  will  be  similar  for  each  of 
the  individual  terms  of  the  above  relationship,  consider 
only  the  specific  term 
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'9y  9y  ' 


where  the  separate  derivatives  can  be  expressed  as 


8ui  _  um  u'(y 1  -  u’(yo) 
3y  y-yc  y  -  yc 


lim 

6y-»-o 


u' 


Combination  of  the  two  quantities  above  then  yields  the 
following : 


, 9u 1  9u 1 
'  9y  9y 


lim 

6y-*-o 


u'u' 

, -+  + 


-  2u'u'  +  u'u* 
- 


The  autocorrelation  portion  can  be  defined  as 


Y  =  7  uju_J_ 

and 

C  =  V u'u' 


By  definition,  the  spacial  correlation  coefficient  f  is 
written  as  (Schlitchting,  1960,  p.  467) : 

uTu1" 

f(T)  =  — - - - 

j  u^_u|  y u '  u ' 

where  f(x)  is  a  function  of  the  separation  distance  t.  I: 
the  definitions  of  equation  5.1  are  used,  the  dissipation 


term  can  be  written  as 


2f  tH'C  + 


dy  dy  x->o 


where  the  separation  distance  t  has  been  substituted  for  6y. 
Two  qualities  of  the  correlation  coefficient  can  now  be 
introduced.  First,  it  can  be  seen  that  as  t  approaches 
zero  f ( t )  approaches  unity.  This  is  merely  a  statement  of 
Cauchy's  inequality  (Taylor,  1955,  p.  204).  Secondly,  it 
is  assumed  that  the  correlation  coefficient  f(x)  is  a 
continuous  function  and,  due  to  Cauchy's  inequality,  is  a 
maximum  at  x  =  0.  By  definition  then,  the  derivative  of 
f(x)  with  respect  to  x,  evaluated  at  the  location  x  =  0, 
must  be  zero.  The  argument  of  a  continuous  function  can  be 
made  by  noting  that  the  quantity  u^_u'  is  finite  when  x  =  0 
and  u+  -+  u.  In  support  of  this  assumption  is  the  often  used 
correlation  coefficient  relationship 


f  (x)  =  e  A 


where  A  is  some  length  parameter  of  the  system.  This  form 
of  relationship  has  been  used  for  isotropic,  homogeneous 
turbulence  as  presented  by  Batchelor  (1956)  and  Hinze  (1959) 
for  a  variety  of  Reynolds  numbers.  The  above  qualities 
can  be  stated  as 


3y  3y 


2f  "cc 


which  reduces  to  the  following  form 


.  3u '  3u' 
l3y  3y 


3  u'u',2 

3y  ; 


f"  u 1 u 1 
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This  relationship  now  expresses  the  dissipation  term 
exactly  in  terms  of  the  Reynolds  stresses.  A  somewhat 
similar  relationship  has  been  developed  for  isotropic  turbu¬ 
lence  (Hinze,  1959,  p.-154)  which  does  not  include  the  first 
term.  It  is  proposed,  therefore,  that  equation  5.5  be  used 
to  approximate  the  autocorrelation  dissipation  terms  in 
the  Reynolds  stress  transport  equations.  The  approximation 
of  the  f"  term  will  be  discussed  at  a  later  point. 

Cross-correlations 

The  cross-correlation  portion  of  the  dissipation  term 
can  be  written  in  the  following  form 

TuTTuT 

As  in  the  previous  development  for  the  autocorrelations, 
consider  just  the  following  term: 

■  3u*  Bv7*. 
l3y  3y  1 
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v*v* ) 


5 


The  above  transformation  is  identical  to  that  shown  by 
Hinze  (1959,  p.  254).  Similar  to  the  autocorrelation 
development,  the  following  definitions  are  presented. 
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—  \J  u*u* 
+  + 

=  J  v*v* 
v  +  + 


A  =  y  u*u* 

b  =  y  v*v* 


and  the  two-point  velocity  correlation  coefficients 


f  ( x )  =  u*u*/  (/ u*u*  /u*u*) 
2_*  +  + 


and 


f  (t)  =  v*v*/  (  y  v*v*  /  v*v* ) 
2*  +  + 


Equation  5.9  can  now  be  written  in  the  following  form  by 
use  of  the  above  definitions: 


9u 1  9v 1 
3y  dy 


lim  1 . 
x->o  2  ' 


i p2  -  62  -  2f.  *\pA  +  2f-*6B  +  A2  -  B 


As  previously  shown,  it  can  be  seen  that  as  x->-o,  ^-*A, 
6-+B,  f^-*-f2->l.  Since  both  the  numerator  and  denominator 
approach  zero  as  t  approaches  zero,  L'Hospital’s  rule  is 
used  to  give 


a 


lizn  1  2'^  "  260  "  2fi*^A  "  2f1*'^’A  +  2f2*6B  +  2f2*6'B 

t+o  2(  2t  * 


since  A  ^  A(t)  and  B  /  B(t).  L'Hospital's  rule  is  again 


used  since  both  numerator  and  denominator  approach  zero 


as  t  approaches  zero.  This  results  in 


(§“  f^-)  =  ^  ^  ( 2ip '  ip '  +  2\pip"  -  2  6*6’  -  266"  -  2f^A 


-  4f.[*(J;,A  -  2f1*^"A  =  2f ^*6B  +  4f^6'B  +  2f2*6"B)/2 


>ince  as  t-*-o;  t|;->A,  6-+B,  f^-fJ^O,  and  f^-^-f^^l,  the  following 


simplification  can  be  made: 


‘If"  If-’  =  I(A'A'  '  B'B’  ‘  £l*AA  +  f2*BB) 


which  reduces  to  the  following  form 


, 9u '  9v ' 

'  9y  9y  1 


1  r  ,9  V  u*u*,  2 

2  1  1  3  t  ’ 


,  9  V  v*v*x  2 
'  9  T  ' 


-  f?*u*u*  +  f?*v*v*] 
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The  above  relationship  can  be  transformed  back  to  the 


original  coordinate  system  using  the  following  relationships: 


.  u’u'  +  2u'v'  +  u'v'  _  .  u'u'  -  2u'v'  +  v'v' 

A)  U*U*  =  - rj -  ,  v*v*  =  - 2 - 


•V-  /  A  A  .SSV  '/'V-  .*• 


•  ***••>  » 


m 


(x*  +  y*) 


C)  since  6t  =  6y  =  6  (- - — Jr-  )  ,  it  can  be  seen  that 
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Substitution  of  these  four  relationships  into  equation  5.10 


gives  the  following 
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The  following  identity  will  be  used  to  obtain  simplified 


forms  of  equations  5.5  and  5.11: 


The  final  forms  for  the  approximation  of  the  dissipation 


of  turbulence  terms  can  be  written  as 


Autocorrelation : 
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Cross-correlation : 
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where  repetition  of  indices  is  not  intended.  The  above 


relationships  are  not  isotropic  and  should  therefore  repre¬ 


sent  an  improvement  over  the  usual  assumption  of  isotropy 


of  dissipation.  This  rigorous  development  was  shown  because 


of  the  relative  importance  of  the  dissipation  terms  in  the 


Reynolds  stress  transport  equations. 
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Pressure  Redistribution 


i  y. 


The  terms  representing  pressure  fluctuation  effects  are 
as  follows 
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where  the  second  group  of  terms  represents  diffusive  trans¬ 
port  and  the  first  group  represents  pressure  redistribution. 
The  diffusion  portion  is  generally  neglected  in  relation  to 
the  redistribution  portion  (as  previously  mentioned, 

Hanjalic  and  Launder,  1972)  in  the  modeling  of  the  terms 
shown  in  equation  5.13.  Results  of  the  order  of  magnitude 
analysis  represented  earlier  confirm  that  other  diffusive 
terms  are  small;  therefore,  in  this  formulation,  only  the 
pressure  redistribution  portion  will  be  considered.  This  is 
the  accepted  practice  in  most  present  turbulence  models. 

An  exact  expression  for  the  first  groups  of  terms  in 
equation  5.13  was  developed  by  Chou  (1945)  based  on  a 
Green's  function  formulation  of  the  Poisson  equation  for 
pressure  fluctuations  (presented  in  Appendix  D) .  Based  on 
Chou's  work,  as  presented  by  Launder,  Reece,  and  Rodi  (1975), 
this  relationship  can  be  written  as 
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where  the  terms  denoted  by  an  asterisk  are  evaluated  at  the 


fixed  vector  point  xQ  and  other  terms  are  evaluated  at  the 


moving  point  xq  +  x.  The  term  represents  a  surface 

integral  which  should  be  negligible  away  from  the  vicinity 


of  the  wall.  The  above  equation  represents  two  separate 


processes,  the  first  involving  just  fluctuating  quantities 


(4>.  .  ,)  and  the  second  involving  the  mean  rate  of  strain 

!j  t  -L 

(<J>.  .  9)  .  Although  equation  5.14  is  an  exact  one,  it  can 
not  be  integrated  and  must  be  approximated.  Unsuccessful 


attempts  to  manipulate  equation  5.14  led  to  the  conclusion 


that  Rotta's  (1951)  proposal  was  the  best  approximation 


available.  It  also  represents  the  most  widely  used  approxi 


mation  for  the  pressure  redistribution  terms. 


Rotta's  proposal  for  the  <J>.  .  0  term  of  equation  5.14 

1  3  r 

can  be  written  as  follows  (as  reported  by  Launder,  Reece, 


and  Rodi ,  1975) : 


’ij,i +  *ji,i  =  -  2/2  6ij  k) 


where  ^  is  a  constant,  and  k  and  e  are  the  turbulent 


kinetic  energy  and  isotropic  dissipation  rate.  According 


to  the  above  reference,  "nearly  every  worker  who  has  made 


closure  approximations  for  (equation  5.14)  has  adopted 


Rotta's  (1951)  proposal  for  4k  ^  In  view  of  this  recom¬ 


mendation,  the  use  of  equation  5.15  seems  fully  justified. 


Rotta's  proposal  for  the  second  term,  cp .  .  was  also 

t  * 


selected  for  use  in  the  present  model.  This  term  can  be 


written  as  follows  (as  presented  by  Launder,  Reece,  and 
Rodi,  1975) : 
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and  £'s  represent  the  Cartesian  components  of  the  position 
vector  x-xq.  According  to  Rotta,  the  fourth-order  tensor 
should  satisfy  the  following  constraints: 


im  im 
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The  form  of  the  tensor  presented  by  Launder,  Reece,  and  Rodi 
(1975)  is  as  follows: 
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where  a,  8,  c 2,  n»  and  v  are  constants.  The  following  rela- 
tionships  result  when  use  is  made  of  the  constraints  of  equa¬ 
tion  5.17. 
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a  =  xy(4c2  +  10) 


6  =  -  yj(2  +  3c2) 


H  =  -  J5(50c2  +  4) 


v  =  '55'(20c2  +  6) 


Although  the  above  relationships  appear  to  be  cumbersome. 


they  are  widely  used  in  the  field  of  turbulence  modeling 


and  have  been  found  to  produce  acceptable  results, 


Diffusion  Closure  Approximation 


The  terms  representing  diffusion  by  turbulent  velocity 


fluctuations  are  the  following  triple  velocity  correlations 


9uiuiuk 


An  exact  transport  equation  for  the  triple  correlation  can 


be  derived  through  the  manipulation  of  equations  C-7  and 


C-9 .  The  full  equation  can  be  written  as 
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As  can  readily  be  seen,  this  equation  contains  many  new 
unknown  correlations.  It  would  be  difficult  to  attempt  to 
scale  these  terms  in  a  manner  similar  to  those  of  the 
previous  equations.  If  it  is  assumed  that  the  time  deriva¬ 
tive,  convective  acceleration,  and  diffusion  terms  are  of 
a  lower  order  of  significance  (as  in  the  Reynolds  stress 
transport  equations) ,  some  simplification  can  be  achieved. 
Further  simplification  can  be  gained  if  use  is  made  of  the 
following  proposal  by  Millionshtchikov  (.1941)  as  reported  by 
Hanjalic  and  Launder  (1972) :  "when  the  triple  correlations 
are  small  and  their  distribution  properties  do  not  differ 
substantially  from  those  of  a  Gaussian  one,  the  quadruple 
correlations  may  be  approximated  in  terms  of  the  second-order 
correlations.  ..."  This  statement  results  in  the  following 
approximation : 


uiujukui  =  uiuj  •  ukui  +  uiuk  •  ujui  +  ului  •  ukuj  5-21 


Use  of  equation  5.21  reduces  the  sum  of  the  quadruple  corre¬ 
lations  and  the  terms  involving  Reynolds  stresses  and  deriva¬ 
tives  of  Reynolds  stresses  to 
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The  above  simplifications  leave  only  the  terms  involving 
pressure  fluctuations  to  be  accounted  for.  If  an  analogy 
exists  between  the  Reynolds  stress  transport  equations  and 
the  equations  of  5.20,  it  cannot  be  assumed  that  the 
pressure  fluctuation  terms  are  insignificant.  In  view  of 
the  above,  and  the  desire  to  use  an  accepted  closure  scheme, 
it  was  decided  to  represent  the  pressure  fluctuation  terms 
by  the  same  approximation  presented  by  Hanjalic  and  Launder 
(1972) .  Their  assumption  of  proportionality  is 
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The  use  of  equations  5.22  and  5.23  gives  the  following 
approximation  for  the  triple  correlation  terms  representing 
diffusion : 
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where  cg  is  a  constant.  Although  severe  simplifications 
were  used  in  the  derivation  of  equation  5.24,  it  is  assumed 
justified  in  view  of  the  order  of  magnitude  analysis  showing 
the  triple  correlations  to  be  of  less  significance  than 
other  terms  in  the  Reynolds  stress  transport  equations. 

The  closure  approximations  presented  in  subsections  5.1, 


5.2,  and  5.3  along  with  the  Reynolds  equations  and  continuity 
equation  now  represent  a  closed  set  of  equations  for 


turbulent  flow.  This  second  order  closure  model  was 
developed  with  a  minimum  of  simplifying  assumptions  so  that 
the  final  set  of  governing  equations  would  be  as  general  as 
possible.  This  requirement  was  felt  to  be  essential  since 
the  ultimate  application  of  the  model  will  be  to  tidally 


induce  flow  regimes. 


SECTION  6 


THE  TURBULENCE  MODEL 


A  consolidation  of  the  previous  sections  can  now  be 
made  so  that  a  completely  closed  set  of  equations  for 
turbulent  flow  of  a  shallow-water  wave  can  be  realized. 

To  accomplish  this  task,  the  Reynolds  equations,  continuity 
equation,  and  Reynolds  stress  transport  equations  must  be 
redimensionalized  and  recombined  to  an  appropriate  and 
equivalent  order  of  magnitude.  The  level  of  accuracy 
selected  for  this  model  was  based  on  the  Reynolds  stress 
transport  equations  such  that  the  lowest  order  of  e  relation 
ship  containing  diffusion  terms  was  included.  This  selec¬ 
tion  was  made  so  that  the  final  set  of  equations  would  be 
parabolic  in  form.  Similar  reasoning  was  made  by  Launder, 
Reece,  and  Rodi  (1975)  in  their  formulation  of  two- 
dimensional  shear  flows.  Inspection  of  equations  4.20 

1/2 

through  4.23  shows  this  level  of  approximation  to  be  0(e  '  ) 
The  equivalent  order  of  accuracy  for  the  Reynolds  equations 
and  the  continuity  equation  is  0(e).  The  terms  including 
these  orders  of  magnitude  will  therefore  represent  the  basic 
turbulence  model. 

In  order  to  explicitly  account  for  a  variable  water- 
surface  boundary,  use  is  made  of  the  continuity  equation 
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and  the  boundary  condition  equations  (presented  in  Appen¬ 
dix  B)  to  express  the  Reynolds  equations  in  a  form  con¬ 
taining  the  water-surface  elevation.  This  is  done  by  first 
recombining  the  Reynolds  equations  and  continuity  equation 
to  an  0(c)  approximation  as  shown  below: 
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Equations  6.1  through  6.3  can  now  be  used  to  develop  the 
governing  equations  in  terms  of  the  water-surface  elevation 
instead  of  the  pressure  term.  The  development  of  this 
begins  with  the  integration  of  the  continuity  equation  (6.3) 
with  respect  to  y  as  shown 


/  IS  +  _/  %  ^  =  0 


where  c  is  the  water-surface  elevation  and  -H  is  the  bottom. 
Integration  can  be  carried  out  by  application  of  Leibnitz's 
rule  for  interchanging  the  order  of  differentiating  and 
integrating  (Taylor,  1955,  p.  523) .  This  yields  the  follow¬ 
ing: 


n mm  mss 


where  the  subscripts  represent  the  mean  velocity  values  at 
the  surface  (s)  and  bottom  (B) .  Use  of  the  boundary  condi¬ 
tion  equations  B.l  and  B.2  reduce  the  above  equation  to  one 
involving  the  water-surface  elevation,  as  shown 
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Approximations  for  the  integral  portion  can  be  made  by 
various  methods,  such  as  the  trapezoidal  rule.  Equation  6.2 
can  also  be  integrated  with  respect  to  y,  for  a  constant 
density  fluid,  as  follows : 


+  g 


dy  + 


9v '  v ' 


9y 


dy 


0 


resulting  in 
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where  -h  is  an  arbitrary  distance  below  the  mean  water  level 
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Since  it  has  been  shown  that  both  — 5 — (C)  and  —5 - (-h)  are 
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0(e  ),  they  are  neglected.  This  relationship  can  be  substi¬ 


tuted  into  equation  6.1  for  the  pressure  term.  The  final 
forms  of  the  Reynolds  equations  and  the  continuity  equation 


to  be  used  in  the  model  can  be  written  in  the  following 
form: 
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Recombination  of  the  Reynolds  stress  transport  equa¬ 
l/2  • 

tions  to  0(e  '  )  results  in  the  following  governing 

equations : 
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The  applicable  closure  approximations,  as  derived  in 
Section  5,  are  as  follows: 
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The  pressure  fluctuation  approximation  terms  can  be  simpli¬ 
fied  by  consideration  of  orders  of  magnitude  of  the  expanded 
terms.  For  instance,  it  has  been  shown  that 
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By  use  of  these  results,  the  pressure  fluctuation  approxi- 
mations  can  be  written  as  follows: 
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where  the  values  of  the  various  coefficients  will  be  dis¬ 


cussed  at  a  later  point. 

Although  this  system  of  equations  appears  to  be  quite 
formidable,  it  does  represent  a  closed  system  of  equations. 
Through  the  judicious  selection  of  an  appropriate  numerical 
solution  scheme  and  an  appropriate  sequence  of  solution, 
equations  6.4-6.10  can  be  solved  for  the  desired  unknown 
quantities  u,  v,  £ ,  u'u',  v'v',  w'w',  and  u'v'.  This  pro¬ 
cedure  will  be  described  in  the  following  section  in 
addition  to  the  applicable  boundary  and  initial  conditions. 
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SECTION  7 


THE  NUMERICAL  APPROXIMATION 

In  this  section,  attention  is  turned  to  the  numerical 
procedure  and  associated  boundary  and  initial  conditions 
necessary  to  arrive  at  a  numerically  stable  approximation 
of  the  governing  equations.  This  must  begin  with  one  of  the 
most  basic  parameters  of  any  numerical  model,  the  vertical 
grid  spacing.  An  equally  spaced  vertical  grid  was  selected 
as  one-tenth  of  the  depth  of  flow.  The  rationale  of  this 
selection  is  as  follows:  inspection  of  the  turbulence  data 
of  Reichardt  (1938)  and  Klebanoff  (1955),  as  presented 
in  Schlichting  (1960,  p.  466,  467) ,  and  the  data  of  Laufer 
(1954) ,  as  presented  in  Hinze  (1959,  p.  521)  show  that  a 
distance  of  one-tenth  of  the  flow  dimension  from  the  wall 
is  outside  the  viscous  sublayer  in  which  the  Reynolds  stress 
values  increase  sharply  from  zero  to  some  maximum  value. 

By  selecting  a  grid  spacing  of  one-tenth  the  flow  depth, 
the  difficulties  of  predicting  this  inflection  point  are 
eliminated.  This  selection  not  only  eases  potential  compu¬ 
tational  difficulties  involving  the  approximation  of  the 
various  terms  but  also  keeps  the  computer  memory  require¬ 
ments  to  a  minimum.  Memory  requirements  and  associated 


computer  run  times  are  an  important  consideration  if  the 
final  model  is  ever  to  be  used  in  a  practical  sense. 

The  numerical  procedure  employed  in  the  model  assumes 
an  initial  mean  velocity  profile  and  solves  for  the  Reynolds 
stresses  for  that  profile  through  an  explicit  iteration 
scheme.  The  mechanics  of  actually  doing  this  involve 
making  initial  estimates  of  the  Reynolds  stresses  and  depth- 
dependent  coefficients  and  updating  those  estimates  during 
each  iteration.  The  resulting  Reynolds  stress  profile  at 
each  x  location  can  then  be  used  to  compute  a  new  mean 
velocity  profile  and  water-surface  profile  for  the  entire 
length  of  flow  at  the  new  time-step  t  +  At.  A  detailed 
step-by-step  description  of  this  procedure  will  be  made 
following  a  description  of  A)  the  initial  and  subsequent 
estimates  of  depth-dependent  coefficients  and  B)  the 
boundary  conditions. 

A)  Depth-Dependent  Coefficients.  Estimates  of  the  depth- 
dependent  parameters  of  equations  6.11-6.22  must  be  made 
before  and  during  the  iteration  procedure.  These  estimates 
are  described  below. 

1.  e/k.  The  value  of  e/k  (or  inversely  k/e)  found  in 
equations  6.11-6.14  and  6.19-6.22  is  computed  from  the 
0(e°)  formulation  of  equations  6.7-6.10  (neglec  ng  the 
diffusion  terms)  as  follows.  By  de£  r  on,  the  isotropic 
dissipation  e  is  expressed  as 
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Since  only  derivatives  with  respect  to  y 
the  pressure  redistribution  terms  vanish 
the  sum  of  equations  6.7-6.10  yields  the 
ship 
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due  to  continuity, 
following  relation- 


u  ’  v’ 


9U 

9y 


+  e  =  0 
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Division  by  the  quantity  k  =  l/2(ukuk)  then  gives  the  value 
of  e/k  for  each  time-step  and  for  each  grid  point  as  shown 


9u 


e/k  =  2  u'v1  —■  /  (u'u'  +  v'v'  +  w'w') 
dy 
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This  coefficient  is  initially  calculated  from  the  input 
Reynolds  stress  and  velocity  profile  data  and  subsequently 
computed  from  the  values  generated  by  the  model.  This 
procedure  avoids  solving  transport  equations  for  c  and  k 


as  is  done  by  Hanjalic  and  Launder  (1972)  and  Launder, 
Reece,  and  Rodi  (1975) . 


2 .  Two-point  velocity  correlation  coefficients.  The  two- 


point  velocity  correlation  coefficient  derivatives  found 
in  the  dissipation  approximations  of  equations  6.15-6.18 
are  also  first  approximated  from  the  0(e°)  formulation  of 
equations  6.7-6.10.  These  first  estimates  are  merely 
calculated  by  balancing  the  sum  of  the  production  and/or 
pressure  redistribution  terms  with  the  diffusion  term. 

This  procedure  was  used  by  Launder,  Reece,  and  Rodi  (1975) 
as  evidenced  by  the  following  excerpt:  "the  initial  value 
of  the  dissipation  rate  was  deduced  as  the  closing  term  in 
a  turbulence  energy  balance."  Subsequent  updatings  of  the 
derivative  of  the  velocity  correlation  coefficient  are  based 
on  the  common  supposition  that  the  dissipation  is  propor¬ 
tional  to  the  turbulent  kinetic  energy  as  shown  by  the 
following  expression  (Rodi,  1980,  p.  21) 


It  is  therefore  assumed  that  each  individual  dissipation 
term  should  be  updated  to  reflect  an  overall  increase  or 
decrease  proportional  to  a  3/2  power  of  the  respective 
Reynolds  stress.  This  is  achieved  by  updating  the  corre¬ 
lation  coefficient  at  the  end  of  each  iteration,  based  on 
the  new  value  of  the  Reynolds  stress  calculated  during 
that  iteration,  using  the  following  relationship: 


/ u ! u ! (t+At) 
f"  (t  +  At)  =  f?  ( T  )  \  /  — — — - — - 

y  u|u|(x) 

where  no  summation  of  indices  is  intended.  This  procedure 
was  found  to  be  necessary  and  effective  in  the  solution  of 
the  Reynolds  stress  transport  equations. 

3.  C*.  The  value  of  C£  in  equation  6.22  was  intended  to 
be  the  constant  in  the  Launder,  Reece,  and  Rodi  (1975) 
formulation.  In  order  to  account  for  wall  effects,  they 
proposed  an  additional  near-wall  correction  term.  Instead 
of  using  their  formulation,  the  following  simpler  correc¬ 
tion  term  was  used: 

C*  =  <c3  +  C3/I) 

where  1  is  the  grid  point  number  measured  from  the  wall, 
i.e.,  1  =  1,  2,  3,  .  .  .,  10.  This  was  found  to  be  quite 
effective. 

4*  */£•  The  value  of  k/e  used  in  equations  6.11-6.13  was 
merely  calculated  as  the  reciprocal  of  e/k  described  in 
section  a)  with  the  following  exception:  because  of  the 
0(e)  formulation  of  e/k  used  in  equation  7.1,  it  can  be 
seen  that  when  either  u 1 v '  or  the  velocity  gradient  approach 
zero,  the  value  of  e  will  approach  zero.  The  reciprocal 
then  approaches  infinity.  This  possibility  sometimes 


occurs  near  the  free  surface  during  the  computation.  Since 
the  term  k/c  only  occurs  in  the  diffusion  term  approxima¬ 
tions  and  it  is  recognized  that  these  are  small  terms, 
especially  near  the  free  surface,  the  following  correction 
was  used: 

C  =  C  (e/k)  when  e/k  <  .40 

s  s 

and 

e/k  =  .20  when  e/k  <  .20 

These  corrections  have  no  theoretical  basis  and  are  only 
used  to  force  the  diffusion  terms  near  zero  close  to  the 
free  surface  where  the  singularity  may  occur.  This  is 
consistent  with  energy  balance  plots  showing  minimal 
diffusion  at  the  surface  (Hinze,  1959) . 

B)  Boundary  Conditions.  Boundary  conditions  used  for  the 
governing  equations  can  be  placed  in  the  following  five 
categories:  (1)  upstream,  (2)  downstream,  (3)  surface, 

(4)  wall,  and  (5)  other.  These  are  addressed  separately 
as  follows. 

1.  Upstream.  The  upstream  boundary  condition  of  known 
discharge  was  maintained;  however,  the  water  level  and 
velocity  profile  were  allowed  to  vary  based  on  the 


following  mass  conservation  considerations.  Equation  6.5 
can  be  written  as 


ii  =  _  la 

3t  3x 
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where  q  is  mass  flux  at  a  given  point  per  unit  mass  and 
width  of  channel.  This  value  is  approximated  by  the 
trapezoidal  rule  as  follows  (Hornbeck,  1975) 


q 


b  -  Av  -  -  m~1- 

f  Udy  =  §2-[U(a)  +  U  (b)  +2  l  U  ( j )  ] 

a  j=l 


where  U(a)  =  0  =  velocity  at  the  wall  and  U(b)  =  surface 
velocity.  Equation  7.3  can  be  written  in  finite  difference 
form  as  follows: 


At  Ax 


Averaging  in  space  and  time  then  gives  the  following 

(q  +  V2  “  -  C°)Ax  =  “(q£  +  q2  -  qi  ~  7.4 


where  superscripts  refer  to  the  new  and  old  time-steps  and 
the  subscripts  refer  to  the  upstream  boundary  and  the  second 
grid  location.  The  new  water  level  can  then  be  calculated 


as  follows: 


q  =  q  ~  q  +  ^2  +  ^qi  +  ql  ”  q2  ~  q^JAt/Ax  7 

A  predictor-corrector  based  new  water  level  was  then 
defined  as 

q  =  l/2(q  +  5°)  7 

Assuming  constant  inflow  discharge,  the  velocity  distribu¬ 
tion  can  be  adjusted  based  on  the  change  in  water  level  as 
follows 

(qq  =  (u°}j  c°/q  7 

where  the  subscript  j  refers  to  the  vertical  grid  points. 

2.  Downstream.  The  downstream  boundary  condition  of  known 
water-level  elevation  was  maintained.  The  velocity  distri 
bution  and  associated  discharge  were  allowed  to  vary 
according  to  an  argument  similar  to  that  used  for  the 
upstream  boundary.  Equation  7.4  can  be  written  for  the 
flux  at  the  N  +  1  (downstream)  boundary  as  follows: 

qN+l  "(;N  +  CN+1  ^N+l)£,X/At  +  qN  +  qN  qN+l 

yielding  the  new  time- step  flux  at  the  downstream  boundary 
The  corresponding  velocity  profile  was  then  calculated 
using  a  similar  profile  assumption  as  follows: 


where  the  subscript  j  refers  to  the  vertical  grid  points. 

3.  Wall.  A  no-slip  boundary  condition  for  the  mean  velocity 
components  and  individual  Reynolds  stresses  was  imposed  at 
the  wall.  In  order  to  calculate  the  derivative  with  respect 
to  y  of  the  Reynolds  stresses  for  the  first  point  from  the 
wall,  however,  a  "wall"  value  was  calculated  for  each 
Reynolds  stress.  These  values  were  then  used  only  for  a 
central  difference  approximation  of  the  derivatives  at  the 
point  adjacent  to  the  wall.  This  "wall"  value  for  the 
normal  stresses  u 1 u' ,  v ' v ' ,  and  w'w'  was  calculated  from  a 
linear  back  extrapolation  of  the  respective  Reynolds  stress 
values  at  the  first  and  second  grid  point  from  the  wall  as 
calculated  from  the  previous  time- step.  This  procedure 
appears  consistent  with  published  data  since  the  first  point 
is  outside  the  viscous  sublayer.  A  further  inspection  of 
this  data  shows  that  the  u’ v'  curve  is  relatively  flat  at 
the  tenth  point,  in  contrast  with  the  normal  Reynolds  strecs 
terms  which  exhibit  a  definite  slope.  For  this  reason,  the 
"wall"  value  for  u ' v '  was  set  equal  to  the  value  calculated 
for  the  first  point  from  the  wall  during  the  previous  time- 
step.  The  derivative  of  u ' v '  was  then  calculated  by  a 
central  difference  approximation  using  the  "wall"  value  and 
the  values  at  the  first  and  second  point  from  the  wall. 


This  condition  also  appears  reasonable  and  provides  accept¬ 
able  results. 


4.  Surface.  The  surface  boundary  condition  imposed  was 
that  of  symmetry  of  the  mean  velocity  profile.  This  con¬ 
dition  forces  the  first  and  second  derivative  of  U  with 
respect  to  y  to  zero  at  the  surface.  An  additional  boundary 
condition  of  u'v'  =  0  was  imposed  at  the  surface.  The 
reasoning  behind  this  is  as  follows:  the  value  of  u'v' 


is  always  less  than  zero  when  the  mean  velocity  gradient  is 


3u 


positive  {■£—  >  0)  .  Conversely,  it  is  greater  than  zero  when 


the  velocity  gradient  is  negative.  This  is  a  consequence 
of  production  being  balanced  partially  by  dissipation,  which 
is  a  positive-definite  quantity.  This  is  verified  by  Hinze 


(1959,  p.  252).  Since  the  sign  of  u ‘ v 1  is  opposite  that 
of  the  velocity  gradient,  and  the  condition  of  symmetry  of 
the  mean  velocity  profile  has  been  imposed,  the  condition 
of  u'v'  =  0  at  the  surface  is  reasonable.  Experimental 
data  also  show  this  to  be  a  reasonable  assumption.  Normal 
Reynolds  stress  values  for  the  surface  were  computed  by 
extrapolation  assuming  a  zero  derivative  with  respect  to 
y  at  the  surface.  This  also  appears  to  be  a  reasonable 
assumption  and  eliminates  the  necessity  of  fixing  the 
Reynolds  stress  values  at  the  surface.  In  addition  to  this 
approximation,  the  parameter  e/k  is  set  equal  to  the  value 
calculated  for  the  point  adjacent  to  the  surface.  This 
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is  of  negligible  consequence  since  the  value  at  that  point 
is  essentially  zero. 


5.  Other.  One  additional  boundary  condition  must  be  imposed 
on  the  Reynolds  stresses.  The  necessity  for  this  can  be 
seen  by  inspection  of  the  Reynolds  stress  transport  equa¬ 
tions  6.7-6.10,  the  associated  approximations  6.11-6.22, 
and  the  imposed  surface  and  bed  boundary  conditions  of 

U  =  V  =  u'u'  =  v ' v '  =  w ' w 1  =  0  at  the  bed  and  3U/9y  = 

2_  2  _ 

3  U/3y  =  u'v'  =  0  at  the  surface.  The  resulting  equations 
then  constitute  a  linear  homogeneous  set  of  simultaneous 
equations,  the  solution  of  which  is  either  the  trivial 
solution  or  an  eigenvalue  solution.  In  order  to  arrive  at 
a  proper  solution,  a  nonzero  boundary  condition  must  be 
imposed.  This  is  determined  from  the  Reynolds  equation  (6.4) 
by  first  writing  the  equation  for  the  surface  as  follows: 


3U 

3t 


+  u  £2 

3x 


+  V 


3U 

3y 


9u  1  v 1 

y  3x  3y 


If  the  condition  of  steady  uniform  flow,  which  is  the 
condition  for  which  this  model  will  be  tested,  is  con¬ 
sidered,  the  resulting  equation  becomes 


3u  1  v 1 
9y 


constant 
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for  a  given  water-surface  slope  corresponding  to  normal 
flow.  Since  the  value  of  u'v'  is  zero  at  the  surface,  the 
additional  boundary  condition  is  calculated,  in  backward 
difference  form,  as  follows 


u 1 v ' k_ i  =  gSAy  7.10 

where  S  is  the  water-surface  slope  and  the  subscript  k-1 
refers  to  the  first  grid  point  below  the  surface. 

Before  a  detailed  sequence-of-events  description  of 
the  model  is  presented,  the  iteration  scheme  used  to 
approximate  the  Reynolds  stresses  is  described.  Each 
Reynolds  stress  equation  (equations  6.7-6.10)  is  first 
expressed  as  the  sum  of  its  respective  production,  dissipa¬ 
tion,  pressure  redistribution,  and  diffusion  segments. 

These  resulting  equations  are  then  written  in  the  form  of 
a  heat  equation  as  expressed  below 


3  u  !  u 
_ i_JL  = 

3t 


- [PROD (u|u j )  +  DISS  (u^u j )  +  REDIS(u|u^)  +  DIFF(uju.!)] 


=  S  UM  ( u  !  u  )  =  0 

l  r 
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where  x  is  some  small  parameter.  A  single  computational 
loop  is  then  written  in  finite  difference  form  containing 
the  following  functional  relationships: 


u'u'(t  +  Ax)  =  u'u'(t)  -  Ax  (SUMuu) 


v'v'(t  +  At)  =  v'v'(t)  -  At (SUMvv) 

7.12 

w'w'(t  +  At)  =  w'w'(t)  -  At (SUMww) 


u'v'(t  +  At)  =  u'v'(t)  =  At (SUMuv) 

The  interation  procedure  then  repeats  this  loop  computing 
a  new  Reynolds  stress  value  at  the  new  t  value  at  the  end 
of  each  iteration.  This  procedure  minimizes  the  sum  terms 
thereby  converging  on  a  stable  Reynolds  stress  value  for 
the  given  boundary  conditions. 

The  following  represents  a  step-by-step  description 
of  the  sequence  of  events  necessary  for  one  completed 
time  iteration  of  the  numerical  model.  It  is  applied  to 
a  section  of  channel  flow  which  is  represented  by  Imax  +  1 
grid  points  in  the  horizontal  direction  and  Kmax  +  1  grid 
points  in  the  vertical.  The  grid  points  corresponding  to 
1=1  represent  the  upstream  boundary,  while  those  corre¬ 
sponding  to  I  =  Imax  +  1  represent  the  downstream  boundary. 
Grid  points  on  the  bottom,  or  wall,  are  respresented  by 
K  =  1,  while  the  surface  is  represented  by  K  =  Kmax  +  1. 
This  represents  Imax  by  Kmax  grid  cells.  The  final  step 
will  then  represent  the  solution  for  the  entire  flow  field 

at  the  time  t  =  t  +  At. 

o 
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Sequence  of  Events: 


1.  The  initial  values  for  the  flow  field  (U^,  £  and  u|uj) 
are  set  equal  to  the  experimental  values. 

2.  The  initial  values  for  the  two-point  velocity  corre¬ 
lation  coefficient  derivatives  are  calculated. 

3.  The  time  loop  begins  for  the  calculation  of  the  flow 

field  for  the  time  t  =  t  +  At  where  t  represents 

o  o 

some  arbitrary  starting  time. 

4.  Mass  flux  at  the  upstream  boundary  is  calculated  for 

-  V 

5.  The  calculation  loop  begins  for  the  horizontal  grid 
points  1=2,  Imax  +  1. 

6.  Initial  values  of  the  parameter  e/k  are  computed  for 
each  vertical  grid  point  K  =  1,  Kmax  +  1  at  the  I  grid 
section. 


Ui 


7.  The  iteration  loop  begins  for  the  calculation  of  the 
Reynolds  stresses  at  the  I  grid  section. 

For  the  dummy  variable  x, 

a)  "wall"  Reynolds  stresses  are  calculated  for  the 
I  grid  section  corresponding  to  K  =  1  for  x. 
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b)  Reynolds  stresses  are  calculated  for  the  remaining 


grid  points  K  =  2,  Kmax  for  x, 

c)  two-point  velocity  correlation  coefficient  deriva¬ 
tives  are  updated  for  each  grid  point  K  =  2, 

Kmax  based  on  the  Reynolds  stress  values  at  t, 

d)  the  parameter  c/k  is  recalculated  using  the  newly 
calculated  values  of  Reynolds  stresses  and  the 
initial  mean  velocity  gradient  for  each  grid  point 
K  =  2,  Kmax  at  x, 

e)  the  normal  Reynolds  stresses  and  the  parameter  e/k 
are  defined  for  the  surface  as  described  under 
boundary  conditions, 

f)  repeat  step  7  for  the  new  iteration  step  x  =  x  +  Ax 

Calculate  the  mean  horizontal  velocity  component  for 
t  =  t  +  At  for  the  grid  points  K  =  2,  Kmax  and  1=2, 
Imax  using  central  difference  approximations. 

Calculate  the  surface  mean  horizontal  velocity  compo¬ 
nent  for  t  =  t  +  At  at  the  grid  point  K  =  Kmax  +  1 
using  backward  difference  approximation  for  the  grid 
sections  1=2,  Imax. 
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10.  Compute  the  water- surf ace  elevation  for  the  grid 

sections  1=2,  Imax  for  t  =  t  +  At. 

o 

11.  Compute  the  mean  velocity  profile  at  the  downstream 
boundary  for  t  =  t  +  At. 

12.  Compute  the  water- surf ace  elevation  for  the  upstream 
boundary  for  t  =  t  +  At. 

13.  Compute  the  mean  velocity  profile  at  the  upstream 
boundary  for  t  =  t  +  At. 

14.  Return  to  step  3  to  begin  calculations  for  the  time- 

step  t  =  t  +  2At. 
o 

These  fourteen  steps  briefly  describe  the  computational 
procedure  involved  in  the  model.  The  next  section  will 
describe  the  actual  application  of  the  model  to  a  set  of 
experimental  data  obtained  in  a  laboratory  flume. 


SECTION  8 


VERIFICATION  OF  THE  MODEL 


Verification  of  the  model  to  experimentally  measured 
data  was  performed  to  demonstrate  that  the  concepts  used  in 
the  derivation  of  the  governing  equations  are  valid.  The 
data  selected  for  this  comparison  are  those  presented  by 
McQuivey  (1973) .  Test  facilities  included  a  20-centimeter 
wide  flume  with  smooth  rigid  boundaries,  both  a  cylindrical¬ 
shaped  and  a  wedge-shaped  hot- film  sensor,  and  a  calibrated 
1/8-inch-diameter  pitot  tube.  Further  details  of  the  equip¬ 
ment  and  measuring  techniques  are  presented  by  McQuivey. 

Mean  flow  parameters  and  variables  for  the  particular  test 
are  shown  in  Tables  1  and  2. 

Table  1 

Mean  Flow  Parameters 


S  =  Slope  of  the  Energy  Gradient 
W  =  Width  of  Flume 
Q  =  Discharge 
T  =  Temperature 

D  =  Depth  of  Flow  (Channel  Center  Line) 
V  =  Average  Velocity  (Q/A) 


=  1.75  x  10  J 
=  0.632  ft. 

=  0.0874  cfs 
=  27 . 2°C 
=  0.1013  ft. 

=  1.36  ft/sec 


R  =  Reynolds  Number  (VD/v) 

F  =  Froude  Number  ( V  >/ g^) 

C/  g  =  Chezy  Discharge  Coefficient 
v  =  Kinematic  Viscosity 


=  15,000 
=  0.75 
=  20.7 

0.920  x  10"5  f t2/sec 


y/D 


U 


u'u' 


v*  v' 


-uV 


>  ^ 


0.775 

1.460 

0.00578 

0.00274 

0.00106 

0.452 

1.350 

0.00922 

0.00417 

0.00251 

0.290 

1.260 

0.01232 

0. 00578 

0.00328 

0.208 

1.190 

0.01513 

0.00711 

0. 00395 

0.129 

1.090 

0.01796 

0.00805 

0.00404 

0.071 

0.960 

0.02560 

0. 01124 

0.00409 

Initial  input  for  the  model  consisted  of  both  the  mean 
flow  variables  for  each  grid  point  and  the  boundary  condi¬ 
tions.  The  mean  flow  values  of  U,  u'u'.  v T v ' ,  and  u' v'  were 
obtained  from  curves  plotted  through  the  experimental  values 
shown  on  Table  2.  The  appropriate  grid  point  value  was  then 
merely  picked  from  the  respective  profile  plot.  Unfortu¬ 
nately,  values  of  w 'w 1  were  not  taken;  therefore,  initial 
estimates  were  made  by  assuming  a  relationship  between 
u'u',  vV,  and  w  'w '  similar  to  that  shown  by  Comte-Bellot '  s 
data  (1965)  as  shown  by  Launder,  Reece,  and  Rodi  (1975)  . 

This  relationship  places  the  w ' w '  value  at  approximately 
one- third  of  the  difference  (u'u'  -  v ' v ' )  above  the  value 
of  v ' v' .  Starting  values  for  w'w'  were  then  estimated  from 
the  plots  of  the  experimental  data.  The  initial  mean  flow 
velocity  profile  and  water-level  elevation  for  the  entire 


flow  section  was  placed  at  a)  V(y)  =  0  (vertical  velocity), 
b)  U(y)  as  plotted  from  Table  2,  and  c)  D(x)  =  0.1013  ft 
(depth  of  flow) .  The  additional  parameter  of  bottom  slope 
had  to  be  deduced  from  the  mean  flow  parameters  (Table  1) 
since  it  was  not  provided  in  the  initial  data.  The  reason 
for  this  is  that  the  slope  of  the  energy  grade  line 
(Table  1)  is  greater  than  the  friction  slope  calculated 
from  the  Chezy  discharge  coefficient,  indicating  gradually 
varying  flow  as  opposed  to  steady  uniform  flow.  A  slightly 
lower  Chezy  discharge  coefficient  of  19.5  was  therefore 
selected  which  yielded  an  intermediate  slope  value  of 
0.00150.  This  value  is  more  realistic  of  the  slope  required 
to  produce  a  normal  depth  of  0.1013  ft,  although  it  is  not 
exact.  This  is  shown  from  the  model  output  which  does 
indicate  a  slight  gradually  varied  flow  condition.  This 
could  be  rectified  by  trial-and-error  determinations; 
however,  it  was  not  felt  necessary  since  the  experimental 
data  also  indicated  gradually  varied  flow. 

Additional  parameters  necessary  for  model  input  are 
those  constants  used  in  the  appropriate  closure  schemes 
mentioned  in  Section  5.  The  resulting  values  were  not 
computer  optimized,  but  were  determined  by  consideration 
of  both  published  values  and  model  output.  Table  3  shows 
a  comparison  of  the  selected  values  with  those  published 


Table  3 


Closure  Approximation  Coefficients 


Symbol 

Present  Model 

Han jalic 
&  Launder 

Launder,  Reece, 
&  Rodi 

C1 

2.00 

2.80 

1.50 

C2 

0.50 

0.  45 

0.40 

C 

s 

0.10 

0.08 

0.11 

C3 

2.50 

“ 

Computational  variables  for  this  selected  flow  regime  are 
as  follows: 

Ax  =  2.0  ft 

Ay  =  1/10  depth  of  flow 
At  =  0.50  sec 
At  =  0.01  sec 

These  values  were  based  on  computational  speed  and  stability 
criteria  and  were  not  computer  optimized.  All  computations 
were  performed  on  a  Texas  Instruments  Model  980-B  mini¬ 
computer  with  64  K  bytes  of  memory.  Peripheral  equipment 
included  a  2.5  megabyte  moving  head  disk,  dual  cassette 
transports,  and  an  electrostatic  printer/plotter.  The 
vertical  mean  velocity  and  Reynolds  stress  distributions 
resulting  from  the  model  are  shown  in  Figures  1  and  2.  The 
mean  velocity  profile  shows  an  inflection  at  the  surface 
which  is  probably  due  to  the  gradually  varied  flow  condition 


Figure  1-  Mean  velocity  profile 


resulting  from  the  channel  bed  slope  selection.  With  this 
exception,  the  agreement  is  acceptable.  The  Reynolds  stress 
profiles  are  in  excellent  agreement  with  the  experimental 
results.  Values  generated  for  w ' w 1  are  included  although 
experimental  results  were  not  provided  for  comparison. 
Figures  3-6  represent  plots  showing  the  computed  distribu¬ 
tion  of  production,  dissipation,  pressure  redistribution, 
and  diffusion  for  each  of  the  Reynolds  stress  equations. 
Figure  7  represents  a  consolidation  of  those  individual 
terms  for  just  the  normal  Reynolds  stress  equations  (Fig¬ 
ures  3-5) .  No  direct  comparisons  can  be  made  to  experi¬ 
mental  data  since  these  terms  are  not  directly  measurable. 

Inspection  of  the  individual  terms  in  the  Reynolds 
stress  transport  equations  (Figures  3-6)  show  some  interest¬ 
ing  trends  which  can  not  be  deduced  from  other  model  results 
For  example,  the  predominant  terms  in  the  equations  for  v'v' 
and  w'w'  are  dissipation  and  pressure  redistribution,  while 
production,  pressure  redistribution,  and  dissipation  are 
significant  in  the  u'u'  equation.  This  would  appear 
reasonable  since  the  v'v'  and  w'w'  equations  contain  no 
production  terms;  therefore  any  turbulent  energy  entering 
from  redistribution  must  be  balanced  by  dissipation. 
Examination  of  the  u'u'  equation  shows  that  the  production 
term  is  balanced  by  a  combination  of  pressure  redistribution 
and  dissipation.  The  u'v'  equation,  which  also  contains  a 
production  type  term,  shows  that  dissipation  is  relatively 
insignificant  and  that  production  is  balanced  by  pressure 


nergy  balance 


-  DISSIPATION 


Figure  7.  energy  balance 
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redistribution.  In  all  cases,  the  diffusion  term  is  small 
as  was  shown  in  the  perturbation  expansions  of  Section  4. 

The  usual  assumption  of  isotropic  dissipation  can  now  be 
evaluated  by  comparing  the  appropriate  terms  in  the  four 
plots.  The  condition  of  negligible  dissipation  in  the  u'v' 
equation  can  be  seen;  however,  dissipation  in  the  u'u' 
equation  is  significantly  lower  than  the  dissipation  in 
the  equations  for  v'v*  and  w'w'.  This  would  lead  to  the 
conclusion,  based  on  this  formulation,  that  turbulent  dissi¬ 
pation  is  not  isotropic  and  that  the  assumption  of  isotropy 
could  lead  to  significant  errors  in  the  overall  model 
results . 

The  above  conclusions  have  pointed  out  that  the 
behavior  of  turbulent  fluctuations  in  a  shallow  water  wave 
climate  is  almost  entirely  a  function  of  the  processes  of 
production,  pressure  redistribution,  and  dissipation.  Of 
these  three  quantities,  only  production  can  be  expressed 
in  terms  of  known  quantities.  The  remaining  terms  must  be 
approximated  through  an  appropriate  closure  scheme.  Sources 
of  error  introduced  through  inaccurate  approximations  are 
difficult  to  detect  due  to  the  high  degree  of  interdepen¬ 
dency  of  the  Reynolds  transport  equations.  Nevertheless, 
an  assessment  of  the  closure  schemes  used  in  this  model 
will  be  made  based  on  the  final  output  of  the  model. 

1.  Diffusion.  Although  it  was  shown  that  the  diffusion 
terms  were  of  a  low  order  of  significance,  they  were 
included  to  render  the  set  of  equations  parabolic.  The 


inclusion  of  these  terms,  especially  near  the  wall,  are 
necessary;  however,  it  is  felt  that  a  simpler  formula¬ 
tion  could  easily  be  used.  The  computational  effort 
required  for  this  approximation  is  much  greater  than 
that  required  for  the  other  closure  schemes.  This  is 
due  not  only  to  the  physical  length  of  .the  approxima¬ 
tion  but  also  to  the  fact  that  the  coefficient  k/e  may 
approach  infinity  near  the  surface.  This  possibility 
must  be  accounted  for.  Due  to  the  relative  magnitude 
of  the  term,  and  the  fact  that  its  accuracy  has  been 
questioned  (Launder,  Reece,  and  Rodi,  1975),  the  above 
observations  lead  to  the  conclusion  that  a  simpler 
approximation  should  be  used. 

Pressure  Redistribution.  The  pressure  redistribution 
terms  are  probably  the  most  important  terms  in  the  set 
of  equations  since  they  are  significant  in  all  equa¬ 
tions.  They  are  also  the  most  difficult  to  approximate 
since  they  deal  with  a  pressure  and  velocity  fluctua¬ 
tion  correlation  term  which  cannot  even  be  measured 
in  the  laboratory.  The  closure  scheme  selected 
(Launder,  Reece,  and  Rodi,  1975)  appears  to  give  very 
reasonable  results.  This  conclusion  is  based  on  both 
the  magnitude  of  the  approximation  and  the  response 
of  the  approximation  to  changes  in  other  flow  vari¬ 
ables.  The  use  of  the  simpler  near-wall  correction 
for  the  u'v'  equation  also  was  shown  to  be  quite 
effective  for  that  region.  In  view  of  the  above,  it 


is  concluded  that  the  pressure  redistribution  closure 
scheme  is  a  very  effective  one.  This  conclusion  was 
also  reached  by  Launder,  Reece,  and  Rodi  (1975). 
Dissipation.  The  dissipation  approximations  presented 
in  this  paper  hopefully  represent  an  improvement  in 
turbulence  modeling  in  that  the  assumption  of  iso¬ 
tropic  dissipation  is  not  made.  The  derivations  for 
the  approximations  are  exact  and  no  major  assumptions 
or  simplifications  were  made.  The  resulting  approxi¬ 
mations  (equations  6.15-6.18)  can  be  simplified  further 
since  the  first  term  is  relatively  insignificant  com¬ 
pared  with  the  term  involving  the  two-point  velocity 
correlation  coefficient  derivatives.  The  normal  dissi¬ 
pation  approximations  then  reduce  to 
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where  summation  of  indices  is  not  intended.  Although 
these  coefficient  derivatives  (f7)  were  initially 
approximate  from  the  experimental  data,  for  this  model 
they  can  be  measured  in  the  laboratory.  This  has  been 
performed  in  a  5-inch  wind  tunnel  by  Laufer  (1951)  using 
the  following  relationship,  identical  with  equation  8.1, 
in  which  he  defines  a  microscale  of  turbulence  for  the 
u'u'  velocity  fluctuation  as 
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where  is  the  appropriate  correlation  coefficient 
for  u ' u ' .  This  is  a  consequence  of  the  assumption  that 


8.3 


Using  the  same  assumption,  equivalent  X  values  can  be 
plotted  for  the  two-point  velocity  correlation  coeffi¬ 
cients  used  in  this  model.  Figure  8  shows  a  comparison 
of  the  model  values  corresponding  to  the  three  Reynolds 
stress  values  in  addition  to  the  data  corresponding  to 
u'u' ,  as  measured  by  Laufer.  Although  absolute  values 
are  not  comparable,  because  of  experimental  differences, 
the  plot  does  demonstrate  that  the  derived  model 
values  are  of  the  same  order  of  magnitude.  This  is 
also  confirmed  by  similar  plots  of  X  values  for  grid 
turbulence  reported  by  Batchelor  and  Townsend  (1948) 
as  shown  in  Cebeci  and  Smith  (1974,  p.  19) .  This  data 
shows  maximum  X  values  of  approximately  0.44.  Although 
the  concept  of  using  velocity  correlation  coefficients 
for  dissipation  approximations  is  not  new,  it  has  not 
been  applied  to  numerical  turbulence  modeling.  In 
view  of  Figure  2,  it  is  apparent  that  it  can  be  success¬ 
fully  used  for  this  purpose  and  does  represent  an 
improvement  over  isotropic  dissipation  models. 
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The  overall  conclusion  regarding  the  verification  of 
the  model  is  that  the  derived  governing  equations  will 
acceptably  reproduce  the  turbulent  flow  regime  for  a  shallow 
water  wave  situation.  This  has  been  shown  through  applica¬ 
tion  of  the  model  to  an  experimentally  derived  set  of  data 
in  which  very  good  reproduction  of  the  laboratory  data  was 
achieved  by  the  model.  The  numerical  methods  employed  in 
this  model  and  the  efficiency  of  these  calculations  were 
not  stressed  since  the  goal  of  this  project  was  to  demon¬ 
strate  the  applicability  of  the  perturbation  approach  for 
a  turbulent  flow  condition  with  practical  applications. 

This  goal  has  been  realized. 


SECTION  9 


CONCLUSIONS 

A  two-dimensional,  second-order  closure,  Reynolds  stress 
turbulence  model  has  been  developed  for  the  investigation  of 
tidally  induced  flow  regimes.  The  governing  equations  for 
this  model  have  been  developed  through  the  application  of  a 
formal  perturbation  expansion  of  the  scaled  and  nondimension- 
alized  turbulent  flow  equations.  Existing  closure  schemes 
were  used  for  the  pressure  redistribution  and  diffusion  terms 
in  the  Reynolds  stress  transport  equations;  however,  a  non¬ 
isotropic  closure  formulation  was  developed  for  the  approxi¬ 
mation  of  the  dissipation  terms.  This  approach  represents  an 
improvement  over  existing  turbulence  models,  which  currently 
use  an  isotropic  dissipation  closure  scheme.  This  feature 
was  developed  in  the  belief  that  turbulent  shear  flows  are 
not  isotropic  and  that  to  consider  them  so  would  severely 
limit  the  versatility  of  the  model.  The  concluding  remark  by 
Launder,  Reece,  and  Rodi  (1975)  that  "...an  assessment  of  the 
validity  of  the  assumption  of  local  isotropy  as  the  most 
urgent  research  task  in  extending  further  the  range  of  appli¬ 
cability  of  the  present  model,"  would  seem  to  reinforce  this 


viewpoint.  This  versatility  is  essential  if  the  model  is  to 
be  applicable  to  complex  flows. 

A  finite  difference  approximation  numerical  model  was 
written  to  solve  the  seven  governing  equations  for  the  mean 
velocities  u  and  v,  the  Reynolds  stresses  u|u 7 ,  and  the 
water-surface  elevation  T .  This  model  was  then  verified  to 
an  existing  set  of  laboratory  data  for  the  case  of  turbulent 
steady  flow  in  a  flume.  The  predicted  mean  and  turbulent 
flow  quantities  agreed  well  with  the  measured  experimental 
values.  Slight  disagreements  may  be  attributed  to  the  fol¬ 
lowing:  1)  the  experimental  data  were  sparse  with  only  six 

measurements  taken  in  the  vertical,  and  2)  steady  uniform 
flow  was  apparently  not  reached  in  the  10-meter-long  flume. 
This  resulted  in  steady,  gradually  varied  flow  in  which  the 
velocity  just  below  the  surface  may  be  greater  than  that  at 
the  surface.  This  was  the  result  shown  by  the  numerical 
model.  In  view  of  the  above,  the  accuracy  of  the  data,  and 
the  numerical  methods  used  in  the  model,  it  is  concluded  that 
the  derived  set  of  governing  equations  and  closure  equations 
do  acceptably  approximate  a  steady  turbulent  flow  regime. 

The  present  investigation  has  shown  the  following: 

1.  The  turbulent  quantity  term  in  the  x-direction  momentum 
equation  (EQ.  6.4)  is  significant  in  magnitude  with 
respect  to  the  mean  flow  quantities.  Improper  modeling 
of  this  term  will  lead  to  inaccuracies  in  the  calcula¬ 
tions  for  mean  flow  quantities. 


Acceleration  terms  and  derivatives  with  respect  to  the 
direction  of  flow  in  the  Reynolds  stress  transport  equa¬ 
tions  are  small  in  comparison  to  other  quantities.  This 
shows  that  changes  of  turbulent  quantities  with  respect 
to  the  vertical  direction  are  more  significant  than 
changes  with  respect  to  the  direction  of  flow.  This 
observation  is  a  result  of  the  fact  that  tidally  induced 
accelerations  in  the  flow  field  are  small.  This  result 
would  not  be  valid  for  steep  waves  causing  a  rapidly 
varying  flow  condition. 

The  nonisotropic  dissipation  closure  approximations 
developed  for  this  project  are  superior,  in  theory,  to 
the  usual  assumption  of  isotropy.  The  approximation 
does  produce  acceptable  results  and  can  be  seen  to  be 
nonisotropic  in  Figures  3-5.  Extensions  of  the  model  to 
nonsteady-state  conditions  should  demonstrate  the  supe¬ 
riority  of  this  closure  approach.  The  potential  impor¬ 
tance  of  this  term  can  also  be  seen  by  the  fact  chat  it 
contains  the  fluid  density  in  the  kinematic  viscosity. 
This  could  be  an  important  consideration  for  future 
extensions  of  the  model  for  density  stratification 
and/or  sediment  transport  studies. 

The  pressure  redistribution  closure  approximation  ap¬ 
pears  to  yield  satisfactory  results.  The  diffusion  clo¬ 
sure  approximation,  however,  has  the  possibility  of 
becoming  unstable  in  areas  of  low  turbulence.  For  this 


reason,  slight  alterations  to  the  present  formulation 
will  probably  have  to  be  made  if  the  model  is  to  be 
tested  for  nonsteady-state  flows. 

The  present  model  was  formulated  to  provide  a  basic 
understanding  of  the  hydrodynamics  of  a  turbulent,  tidally 
induced  flow.  Without  this  insight,  it  would  be  difficult  to 
conduct  a  meaningful  investigation  of  other  flow-related 
phenomena  such  as  density  stratification  or  the  nonlinear 
propagation  of  tides.  The  model  developed  for  this  project 
has  been  shown  to  adequately  represent  the  various  mechanisms 
of  steady  turbulent  flow.  The  extension  to  nonsteady  flow 
should  not  be  too  difficult  since  the  basic  formulation  was 
for  that  specific  application.  Following  verification  to 


nonsteady  conditions,  the  model  can  be  used  in  conjunction 
with  additional  equations  for  the  detailed  study  of  other 
flow-related  problems.  An  additional  benefit  of  the  model  is 
that  the  concepts  developed  in  the  formulation  of  the  model 
can  be  used  to  improve  upon  existing  hydrodynamic  models.  In 
either  application,  the  model  will  have  produced  an  increased 
understanding  of  turbulent  flow  which  will  ultimately  aid  in 
the  development  of  improved  predictive  models  for  estuarine 
flow. 
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APPENDIX  A 


THE  REYNOLDS  EQUATIONS  AND  THE  CONTINUITY  EQUATION 

The  equations  governing  the  flow  distribution  of  a 
fluid  are  based  on  the  principle  of  conservation  of  mass 
and  conservation  of  linear  momentum.  The  two-dimensional 
equations  for  an  incompressible,  Newtonial  fluid  can  be 
written  as  follows: 


Conservation  of  momentum:  Navier-Stokes  equations 


H  +  uBI  + 

at  ax  ay 


i  w  yu  aV 

"  8*  ax2  ay2’ 


A.  1 


9V  +  ulV  3V 
3t  3x  v3y 


1  3P  .  ,32V  .  32V. 


A. 2 


Conservation  of  mass:  continuity  equation 


3U  3V 
3x  3y 


0 


A. 3 


where  U  and  V  represent  the  instantaneous  horizontal  and 
vertical  velocities,  P  represents  the  instantaneous  pressure 
g  is  the  negative  body  force  due  to  gravity,  and  p  and  v 
represent  the  fluid  density  and  kinematic  viscosity.  The 


above  equations  correctly  describe  turbulent  motion;  however, 
the  time  and  length  scales  necessary  to  do  so  are  beyond 
the  practical  capabilities  of  present  computers.  The 
accepted  approach  to  describing  the  turbulent  motion  was 
first  proposed  by  Osborn  Reynolds.  The  suggested  approach 
was  that  each  flow  variable  should  be  represented  as  the 
sum  of  an  average  quantity  plus  a  fluctuating  quantity. 

This  relationship  can  be  expressed  as  follows 

U  =  U  +  u' 

V  «  V  +  v'  A. 4 

P  =  P  +  p' 

where  the  average  quantities  are  defined  as 

T  T  T 

U  =  i  /  U  (t+T )  dx ,  V  =  |  /  V(t+T)dT,  P  =  i  f  P(t+T)dT 

o  o  o 

The  magnitude  of  T  is  further  defined  as  large  in  comparison 
to  the  time  scales  of  turbulent  motion.  Substitution  of 
equations  A. 4  into  equations  A.l,  A. 2,  and  A. 3,  and  averaging 
the  resulting  quantities  with  respect  to  time  yields  the 
following  equations ; 


m 


1Y  +  g3V  +  -9V  = 
3t  3x  v  3y 


1  3P  .  ,  32V  ^  32V. 

O  ?77  -  g  +  - o  +  y) 

p  3y  3x2  3y2 


3u* v1  _  9v' v1 
3x  3y 


+  3V  = 

3x  3y  u 


Equations  A. 5  and  A. 6  represent  the  well-known  Reynolds 
equations  and  equation  A. 7  is  the  continuity  equation.  The 
fluctuating  quantities  are  referred  to  as  the  Reynolds 
stresses.  The  goal  of  a  turbulent  model  is  to  provide  a 
solution  to  these  equations  for  the  mean  flow  quantities. 
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APPENDIX  b 
boundary  conditions 


•  .  £or  the  free  surface  and  bottom  of 

— ;;n;:::cmed  based  on  the  dimensions  shown 
the  flow  region  are  p 
below  (for  two-dimensional  flow) 


f  sur£aoe  u  defined  as  a  wave  superimposed  on 

where  the  free  developed  by 

,er  level.  Boundary  conditions  are 

*  me“  .  bl.  F(x  y.t,  for  both  the  upper  and  lower 

defining  a  variable  Flx.y.ti 

.  has  the  following  property 

surface  which  has 


dF  -  0  on  the  boundary 
dt  " 


The  surface  variable  can  be  defined  as 


P  m  y  -  Ux,t) 

s 


Bl 


The  boundary  equation  can  then  be  written  as  follows 


»'LS: 
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dF  3F  3F„  3F 

dt“  =  lt~  +  Us  Jt~  +  Vs  3y“  =  0 


If  +  Us  If  +  Vs  -  0 


where  U  and  V  represent  velocities  in  the  x  and  y  directions 
and  the  subscript  denotes  the  surface.  The  bottom  variable 
can  be  defined  as  shown  below 


Fb  »  y  +  H(x) 


The  resulting  equation  can  be  written  as 


A  turbulent  expansion,  similar  to  that  shown  in  Appendix  A, 
yields  the  following  equations 


.iiL  +  u  +  u '  -  v  =  0 

3t  us  3x  us  3x  vs 


u  iii  +  v 

UB  3x  VB 


These  two  relationships  define  the  flow  boundaries  for 
turbulent  free  surface  flow  in  two  dimensions. 
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APPENDIX  C 


THE  REYNOLDS  STRESS  TRANSPORT  EQUATIONS 


The  transport  equations  for  the  Reynolds  stresses 
(uTuT)  can  be  derived  from  the  Navier-S tokes  equations 
shown  in  equations  A.l  and  A. 2.  This  can  be  accomplished 
by  first  writing  the  equation  in  indicial  notation  as 
follows : 


3U .  3U.  .  .  32U. 

_ i.  ,  ri  i _ 1  9P  .  «  +  ^ 

3 t  uk  3xk  “  p  3x±  ri  9xk3xk 


C.l 


where  F.  represents  the  body  force  (F  =  0,  F  =  -  g) .  As 
x  x  y 

was  shown  in  Appendix  A,  each  variable  is  expressed  as  the 
sum  of  an  average  quantity  plus  a  fluctuating  quantity  as 
shown  below 


P  =  P  +  p» 

Substitution  of  the  expansions  of  equations  C.2  into 
equation  C.l  gives  the  following 


Cl 


rig 


X,  dX, 

k  k 


P  dXi 


dxk3xk 


Substitution  of  equation  C.5  from  C.4  results  in  an  equation 
for  the  transport  of  the  single  velocity  fluctuation  u| . 

This  relationship  becomes: 


k  3x^  3x^  1  k  i  k  p  3x^  dx^dx^ 


This  can  also  be  written  for  the  j  index  as: 


+  uk  33^  +  7^  (ujuk  '  ujuk  >  = 


-  I|El+  v 

P  SXj 


multiplication  of  equation  C.6  by  u^ ,  multiplication  of 
equation  C.7  by  u|,  and  addition  of  the  two  results  become 

3u!  3u’.  3U .  3U .  _  3u!  _  3u 

uj  IT-  +  ui  TE1  +  ujuk  +  uiuk  +  ujUk  ^  +  uiUk  9x 


+  uj  3^  (uiuk  -  uiui  >  +  ui  3^  tujuk  -  ujuk  > 


32u! 


3  ui 


“  ■  P  uj  Hr  "  Q  ui  15“  +  v  uj  3xk3xk  *  v  ui  3xkLk 


The  following  identities  can  be  seen: 


ulU.  ■5-^=.  +  u '  U,  ^r-^-  =  u,  - 

3  k  3x^  1  k  3x^  k  3xk 


Substitution  of  these  relationships  into  equation  C.8  and 
averaging  the  results  yields  the  following: 


3 u !  u _  3U.  _  3U.  _  3u u - r - 

+  u  '■  u,'  ^  +  u  !  U '  -s— i  +  U,  ■-1-3  +  u TT^—  u !  u,' 

3t  j  k  3x,  i  k  3x,  k  3x,  i  3x,  l  k 


uu  .  u, 
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1  9xk 
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P  U 

i  3x_.  uj  3x,3x 

j  x 

32“5 

+  V  u!  -5 - ± 

1  3xk 

The  following 

expansions 

can 

now  be  introduced: 

3  u  :  u u 


+  uiuk  "3x 


3  u:  u 


i  3x,  3x 


Substitution  of  these  expansions  into  equation  C.9  results 
in  the  following  familiar  form  of  the  Reynolds  stress 
transport  equations: 
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II 
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p-  8u1  p.  1  3p'u!  1  3p'u! 

p  3x.  —  p  3x  .  p  3x.  p  3x . 


III 


IV 


+  v 
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3u!  3  u!.  Su'.ulu' 

1  1  1  ]  K 

3x,  3x,  3x, 

k  k  k 

VI  VII 


C.10 


The  usual  interpretation  of  the  various  terms  is  presented 
below: 

I :  These  two  terms  represent  the  total  rate  of  change  of 

the  turbulent  fluctuation  (d/dt) .  The  first  represents 
the  time  rate  of  change  while  the  second  represents  the 
convection  of  the  turbulent  fluctuation. 


II:  These  terms,  usually  referred  to  as  "production"  terms, 

represent  the  transfer  of  turbulent  energy  in  the 
system  by  providing  a  link  between  mean  flow  quantities 
and  turbulent  quantities  (Reynolds  stresses) . 

Ill:  Terms  III  are  usually  referred  to  as  the  "pressure 

redistribution"  terms.  This  pressure  rate-of-strain 
correlation  has  the  effect  of  transferring  energy  in 


C5 


the  system  to  a  more  isotropic  condition.  When  complete 
isotropy  is  achieved,  it  can  be  seen  that  this  term 
vanishes  due  to  the  continuity  equation  (i.e.,  when 

i  =  j)  • 

IV:  Terms  IV  are  referred  to  as  pressure-diffusion  terms. 

They  are  generally  much  smaller  in  magnitude  than  the 
pressure  redistribution  terms  and  are  therefore  usually 
neglected. 

V:  These  terms  represent  diffusion  of  the  Reynolds  stresses 

through  viscous  action.  They  are  also  generally  small 
in  magnitude. 

VI:  The  terms  labeled  VI  are  referred  to  as  the  "dissipa¬ 

tion"  terms.  They  are  responsible  for  dissipating 
turbulent  energy  through  the  viscous  interactions  of 
the  rate  of  strain  correlations. 

VII:  The  triple  correlations  of  term  VII  represent  the 
convection  of  the  turbulent  fluctuations  by  the 
turbulent  fluctuations.  This  term  will  be  referred  to 
as  the  "diffusion"  term. 

Equations  C.10  are  the  governing  equations  for  the 
transport  of  the  turbulent  Reynolds  stresses.  The  solution, 
or  approximation,  of  these  equations  is  the  ultimate  goal  of 
a  turbulence  model. 


APPENDIX  D 


THE  POISSON  EQUATION  FOR  PRESSURE  FLUCTUATIONS 

The  divergence  of  the  Navier-Stokes  equations  yields 
a  Poisson  equation  for  the  pressure  term  for  the  case  of 
incompressible  flow.  This  relationship  can  be  written  for 
the  two-dimensional  case  as  follows: 

1  02  P  _  32UU  ,  „  32UV  ,  32W 
-  p  7  P  =  3x2  3*3y  3y2 

As  was  shown  in  Appendix  A,  each  variable  can  be  expressed 
as  the  sum  of  an  average  quantity  plus  a  turbulent  fluctuat¬ 
ing  quantity,  as  shown  below 

U  =  U  +  u' 

V  =  V  +  v' 

P  =  P  +  p' 


Substitution  of  these  expansions  into  equation  D.l  yields 
the  following: 


a  u'u 


9  u  v 


a  v  v 


If  equation  D.l  is  averaged  and  then  subtracted  from 


equation  D.2,  the  resulting  equation  is  a  Poisson  equation 


for  the  pressure  fluctuation.  This  relationship  can  be 


written  as  follows 


2  2 

I  v2  p'  =  (2Uu'  +  u'u')  +  2  I-*-  (Uv'  +  Vu'  +  u'v') 

P  dXdy 


+  -2-^5-  (2Vv‘  +  v '  v '  ) 

3y 


Equation  D. 3  can  be  rearranged,  by  use  of  the  continuity 


equation,  and  expressed  in  the  following  indicial  notation: 


92u!u'.  3U .  9u'. 

_ — i _ X  +  2  i _ i-  _ 2- 

9x^9Xj  9Xj  9x^ 


p  3xk9xk 


The  above  equation  is  that  used  by  Chou  (1945)  for  the  exact 


form  of  the  pressure  fluctuation  -  rate  of  strain  correla¬ 


tion  shown  in  equation  5.14.  Although  a  closed- form  solu¬ 


tion  is  not  possible,  the  use  of  equation  D.4  has  led  to 


satisfactory  approximations  for  the  pressure  fluctuation 


term  in  the  Reynolds  stress  transport  equations. 
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